import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
from torch.utils.data import Dataset
from torch.utils.data import DataLoader
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
from time import sleep
import os
from os import listdir
import matplotlib.image as mpimg
from matplotlib.animation import FuncAnimation as FA
import random
pd.set_option('display.max_columns', 50)
pd.set_option('display.max_rows', 50)
Setting seeds for ensuring reproducibility
# import time
# # Generate a random seed based on the current time
# seed = int(time.time())
# print("Selected seed:", seed)
seed = 42
torch.manual_seed(seed)
torch.cuda.manual_seed_all(seed)
np.random.seed(seed)
random.seed(seed)
📌 Loading the data:
There are 100 trajectories of engine degradation in EACH of the datasets (train and test).
For the train data, each engine is run from a certain normal condition till failure.
For the test data, each engine is NOT run until failure (or we have data until a specific point of the engine's state). The amount of time cycles left for this engine to still be 'normal' is what we call RUL, or Remaining Useful Lifetime.
Predicting what the RUL is for the last state of each machine in the test set will be our prediction task.
folder_path = './CMAPSSData/'
listdir(folder_path)
file_name = 'FD004.txt'
df_train = pd.read_csv(folder_path + 'train_' + file_name, header = None, sep = ' ')
df_test = pd.read_csv(folder_path + 'test_'+file_name, header = None, sep = ' ')
rul_test = pd.read_csv(folder_path + 'RUL_'+file_name, header = None)
for df in [df_train, df_test, rul_test]:
display(df.head())
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 42.0049 | 0.8400 | 100.0 | 445.00 | 549.68 | 1343.43 | 1112.93 | 3.91 | 5.70 | 137.36 | 2211.86 | 8311.32 | 1.01 | 41.69 | 129.78 | 2387.99 | 8074.83 | 9.3335 | 0.02 | 330 | 2212 | 100.00 | 10.62 | 6.3670 | NaN | NaN |
| 1 | 1 | 2 | 20.0020 | 0.7002 | 100.0 | 491.19 | 606.07 | 1477.61 | 1237.50 | 9.35 | 13.61 | 332.10 | 2323.66 | 8713.60 | 1.07 | 43.94 | 312.59 | 2387.73 | 8046.13 | 9.1913 | 0.02 | 361 | 2324 | 100.00 | 24.37 | 14.6552 | NaN | NaN |
| 2 | 1 | 3 | 42.0038 | 0.8409 | 100.0 | 445.00 | 548.95 | 1343.12 | 1117.05 | 3.91 | 5.69 | 138.18 | 2211.92 | 8306.69 | 1.01 | 41.66 | 129.62 | 2387.97 | 8066.62 | 9.4007 | 0.02 | 329 | 2212 | 100.00 | 10.48 | 6.4213 | NaN | NaN |
| 3 | 1 | 4 | 42.0000 | 0.8400 | 100.0 | 445.00 | 548.70 | 1341.24 | 1118.03 | 3.91 | 5.70 | 137.98 | 2211.88 | 8312.35 | 1.02 | 41.68 | 129.80 | 2388.02 | 8076.05 | 9.3369 | 0.02 | 328 | 2212 | 100.00 | 10.54 | 6.4176 | NaN | NaN |
| 4 | 1 | 5 | 25.0063 | 0.6207 | 60.0 | 462.54 | 536.10 | 1255.23 | 1033.59 | 7.05 | 9.00 | 174.82 | 1915.22 | 7994.94 | 0.93 | 36.48 | 164.11 | 2028.08 | 7865.80 | 10.8366 | 0.02 | 305 | 1915 | 84.93 | 14.03 | 8.6754 | NaN | NaN |
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 20.0072 | 0.7000 | 100.0 | 491.19 | 606.67 | 1481.04 | 1227.81 | 9.35 | 13.60 | 332.52 | 2323.67 | 8704.98 | 1.07 | 43.83 | 313.03 | 2387.78 | 8048.98 | 9.2229 | 0.02 | 362 | 2324 | 100.00 | 24.31 | 14.7007 | NaN | NaN |
| 1 | 1 | 2 | 24.9984 | 0.6200 | 60.0 | 462.54 | 536.22 | 1256.17 | 1031.48 | 7.05 | 9.00 | 174.46 | 1915.21 | 7999.94 | 0.93 | 36.11 | 163.61 | 2028.09 | 7863.46 | 10.8632 | 0.02 | 306 | 1915 | 84.93 | 14.36 | 8.5748 | NaN | NaN |
| 2 | 1 | 3 | 42.0000 | 0.8420 | 100.0 | 445.00 | 549.23 | 1340.13 | 1105.88 | 3.91 | 5.69 | 137.34 | 2211.93 | 8305.38 | 1.01 | 41.52 | 129.98 | 2387.95 | 8071.13 | 9.3960 | 0.02 | 328 | 2212 | 100.00 | 10.39 | 6.4365 | NaN | NaN |
| 3 | 1 | 4 | 42.0035 | 0.8402 | 100.0 | 445.00 | 549.19 | 1339.70 | 1107.26 | 3.91 | 5.69 | 137.23 | 2211.89 | 8301.00 | 1.01 | 41.73 | 129.48 | 2387.90 | 8078.89 | 9.3594 | 0.02 | 328 | 2212 | 100.00 | 10.56 | 6.2367 | NaN | NaN |
| 4 | 1 | 5 | 35.0079 | 0.8400 | 100.0 | 449.44 | 555.10 | 1353.04 | 1117.80 | 5.48 | 7.97 | 192.94 | 2222.71 | 8331.05 | 1.02 | 41.32 | 181.82 | 2387.87 | 8057.83 | 9.3030 | 0.02 | 333 | 2223 | 100.00 | 14.85 | 8.9326 | NaN | NaN |
| 0 | |
|---|---|
| 0 | 22 |
| 1 | 39 |
| 2 | 107 |
| 3 | 75 |
| 4 | 149 |
📌 Attaching column names: We have three operational setting columns (os + number), and 21 sensor columns (s + number). We have dropped the last two columns containing NaNs.
col_names = []
col_names.append('unit')
col_names.append('time')
for i in range(1,4):
col_names.append('os'+str(i))
for i in range(1,22):
col_names.append('s'+str(i))
df_train = df_train.iloc[:,:-2].copy()
df_train.columns = col_names
display(df_train.head())
df_test = df_test.iloc[:,:-2].copy()
df_test.columns = col_names
display(df_test.head())
| unit | time | os1 | os2 | os3 | s1 | s2 | s3 | s4 | s5 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s18 | s19 | s20 | s21 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 42.0049 | 0.8400 | 100.0 | 445.00 | 549.68 | 1343.43 | 1112.93 | 3.91 | 5.70 | 137.36 | 2211.86 | 8311.32 | 1.01 | 41.69 | 129.78 | 2387.99 | 8074.83 | 9.3335 | 0.02 | 330 | 2212 | 100.00 | 10.62 | 6.3670 |
| 1 | 1 | 2 | 20.0020 | 0.7002 | 100.0 | 491.19 | 606.07 | 1477.61 | 1237.50 | 9.35 | 13.61 | 332.10 | 2323.66 | 8713.60 | 1.07 | 43.94 | 312.59 | 2387.73 | 8046.13 | 9.1913 | 0.02 | 361 | 2324 | 100.00 | 24.37 | 14.6552 |
| 2 | 1 | 3 | 42.0038 | 0.8409 | 100.0 | 445.00 | 548.95 | 1343.12 | 1117.05 | 3.91 | 5.69 | 138.18 | 2211.92 | 8306.69 | 1.01 | 41.66 | 129.62 | 2387.97 | 8066.62 | 9.4007 | 0.02 | 329 | 2212 | 100.00 | 10.48 | 6.4213 |
| 3 | 1 | 4 | 42.0000 | 0.8400 | 100.0 | 445.00 | 548.70 | 1341.24 | 1118.03 | 3.91 | 5.70 | 137.98 | 2211.88 | 8312.35 | 1.02 | 41.68 | 129.80 | 2388.02 | 8076.05 | 9.3369 | 0.02 | 328 | 2212 | 100.00 | 10.54 | 6.4176 |
| 4 | 1 | 5 | 25.0063 | 0.6207 | 60.0 | 462.54 | 536.10 | 1255.23 | 1033.59 | 7.05 | 9.00 | 174.82 | 1915.22 | 7994.94 | 0.93 | 36.48 | 164.11 | 2028.08 | 7865.80 | 10.8366 | 0.02 | 305 | 1915 | 84.93 | 14.03 | 8.6754 |
| unit | time | os1 | os2 | os3 | s1 | s2 | s3 | s4 | s5 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s18 | s19 | s20 | s21 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 20.0072 | 0.7000 | 100.0 | 491.19 | 606.67 | 1481.04 | 1227.81 | 9.35 | 13.60 | 332.52 | 2323.67 | 8704.98 | 1.07 | 43.83 | 313.03 | 2387.78 | 8048.98 | 9.2229 | 0.02 | 362 | 2324 | 100.00 | 24.31 | 14.7007 |
| 1 | 1 | 2 | 24.9984 | 0.6200 | 60.0 | 462.54 | 536.22 | 1256.17 | 1031.48 | 7.05 | 9.00 | 174.46 | 1915.21 | 7999.94 | 0.93 | 36.11 | 163.61 | 2028.09 | 7863.46 | 10.8632 | 0.02 | 306 | 1915 | 84.93 | 14.36 | 8.5748 |
| 2 | 1 | 3 | 42.0000 | 0.8420 | 100.0 | 445.00 | 549.23 | 1340.13 | 1105.88 | 3.91 | 5.69 | 137.34 | 2211.93 | 8305.38 | 1.01 | 41.52 | 129.98 | 2387.95 | 8071.13 | 9.3960 | 0.02 | 328 | 2212 | 100.00 | 10.39 | 6.4365 |
| 3 | 1 | 4 | 42.0035 | 0.8402 | 100.0 | 445.00 | 549.19 | 1339.70 | 1107.26 | 3.91 | 5.69 | 137.23 | 2211.89 | 8301.00 | 1.01 | 41.73 | 129.48 | 2387.90 | 8078.89 | 9.3594 | 0.02 | 328 | 2212 | 100.00 | 10.56 | 6.2367 |
| 4 | 1 | 5 | 35.0079 | 0.8400 | 100.0 | 449.44 | 555.10 | 1353.04 | 1117.80 | 5.48 | 7.97 | 192.94 | 2222.71 | 8331.05 | 1.02 | 41.32 | 181.82 | 2387.87 | 8057.83 | 9.3030 | 0.02 | 333 | 2223 | 100.00 | 14.85 | 8.9326 |
📌 Attaching RUL(remaining useful lifetime) values to the datasets.
For the train data, the RUL values are not specified, but the document regarding this dataset specifies that for the training data, all the engines were run to failure. Thus, for example, if we had five rows for a specific unit, say unit 7:
| unit | time cycle |
|---|---|
| 7 | 1 |
| 7 | 2 |
| 7 | 3 |
| 7 | 4 |
| 7 | 5 |
Then we know that the last row is when the RUL value becomes 0 (failure), so the RUL for this unit would be attached in this way:
| unit | time cycle | RUL |
|---|---|---|
| 7 | 1 | 4 |
| 7 | 2 | 3 |
| 7 | 3 | 2 |
| 7 | 4 | 1 |
| 7 | 5 | 0 |
For the test data, we have the 'solutions' for the test engines in a separate file called rul_test (the name that I used for dataframe). I will use that to attach the RUL values for the test as well.
units_training = max(df_train['unit'])
units_testing = max(df_test['unit'])
print(f'Units in the training dataset: {units_training}')
print(f'Units in the testing dataset: {units_testing}')
df_train = df_train[df_train['unit'] != 260]
if units_training > units_testing:
df_train = df_train[df_train['unit'] <= units_testing]
elif units_training < units_testing:
df_test = df_test[df_test['unit'] <= units_training]
assert max(df_train['unit']) == max(df_test['unit'])
no_units = max(df_train['unit'])
print(f'No of units in training and testing dataset after: {no_units}')
Units in the training dataset: 249 Units in the testing dataset: 248 No of units in training and testing dataset after: 248
MAX_RUL = 125
no_units = min(rul_test.shape[0],max(df_train['unit']))
print(f'units :{no_units}')
units :248
rul_list = []
engine_numbers = no_units
for n in np.arange(1,engine_numbers+1):
time_list = np.array(df_train[df_train['unit'] == n]['time'])
length = len(time_list)
rul = list(length - time_list)
rul = [min(MAX_RUL,x) for x in rul]
rul_list += rul
df_train['rul'] = rul_list
rul_list = []
for n in np.arange(1,engine_numbers+1):
time_list = np.array(df_test[df_test['unit'] == n]['time'])
length = len(time_list)
rul_val = rul_test.iloc[n-1].item()
rul = list(length - time_list + rul_val)
rul = [min(MAX_RUL,x) for x in rul]
rul_list += rul
df_test['rul'] = rul_list
for df in [df_train, df_test]:
display(df.head())
| unit | time | os1 | os2 | os3 | s1 | s2 | s3 | s4 | s5 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s18 | s19 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 42.0049 | 0.8400 | 100.0 | 445.00 | 549.68 | 1343.43 | 1112.93 | 3.91 | 5.70 | 137.36 | 2211.86 | 8311.32 | 1.01 | 41.69 | 129.78 | 2387.99 | 8074.83 | 9.3335 | 0.02 | 330 | 2212 | 100.00 | 10.62 | 6.3670 | 125 |
| 1 | 1 | 2 | 20.0020 | 0.7002 | 100.0 | 491.19 | 606.07 | 1477.61 | 1237.50 | 9.35 | 13.61 | 332.10 | 2323.66 | 8713.60 | 1.07 | 43.94 | 312.59 | 2387.73 | 8046.13 | 9.1913 | 0.02 | 361 | 2324 | 100.00 | 24.37 | 14.6552 | 125 |
| 2 | 1 | 3 | 42.0038 | 0.8409 | 100.0 | 445.00 | 548.95 | 1343.12 | 1117.05 | 3.91 | 5.69 | 138.18 | 2211.92 | 8306.69 | 1.01 | 41.66 | 129.62 | 2387.97 | 8066.62 | 9.4007 | 0.02 | 329 | 2212 | 100.00 | 10.48 | 6.4213 | 125 |
| 3 | 1 | 4 | 42.0000 | 0.8400 | 100.0 | 445.00 | 548.70 | 1341.24 | 1118.03 | 3.91 | 5.70 | 137.98 | 2211.88 | 8312.35 | 1.02 | 41.68 | 129.80 | 2388.02 | 8076.05 | 9.3369 | 0.02 | 328 | 2212 | 100.00 | 10.54 | 6.4176 | 125 |
| 4 | 1 | 5 | 25.0063 | 0.6207 | 60.0 | 462.54 | 536.10 | 1255.23 | 1033.59 | 7.05 | 9.00 | 174.82 | 1915.22 | 7994.94 | 0.93 | 36.48 | 164.11 | 2028.08 | 7865.80 | 10.8366 | 0.02 | 305 | 1915 | 84.93 | 14.03 | 8.6754 | 125 |
| unit | time | os1 | os2 | os3 | s1 | s2 | s3 | s4 | s5 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s18 | s19 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 20.0072 | 0.7000 | 100.0 | 491.19 | 606.67 | 1481.04 | 1227.81 | 9.35 | 13.60 | 332.52 | 2323.67 | 8704.98 | 1.07 | 43.83 | 313.03 | 2387.78 | 8048.98 | 9.2229 | 0.02 | 362 | 2324 | 100.00 | 24.31 | 14.7007 | 125 |
| 1 | 1 | 2 | 24.9984 | 0.6200 | 60.0 | 462.54 | 536.22 | 1256.17 | 1031.48 | 7.05 | 9.00 | 174.46 | 1915.21 | 7999.94 | 0.93 | 36.11 | 163.61 | 2028.09 | 7863.46 | 10.8632 | 0.02 | 306 | 1915 | 84.93 | 14.36 | 8.5748 | 125 |
| 2 | 1 | 3 | 42.0000 | 0.8420 | 100.0 | 445.00 | 549.23 | 1340.13 | 1105.88 | 3.91 | 5.69 | 137.34 | 2211.93 | 8305.38 | 1.01 | 41.52 | 129.98 | 2387.95 | 8071.13 | 9.3960 | 0.02 | 328 | 2212 | 100.00 | 10.39 | 6.4365 | 125 |
| 3 | 1 | 4 | 42.0035 | 0.8402 | 100.0 | 445.00 | 549.19 | 1339.70 | 1107.26 | 3.91 | 5.69 | 137.23 | 2211.89 | 8301.00 | 1.01 | 41.73 | 129.48 | 2387.90 | 8078.89 | 9.3594 | 0.02 | 328 | 2212 | 100.00 | 10.56 | 6.2367 | 125 |
| 4 | 1 | 5 | 35.0079 | 0.8400 | 100.0 | 449.44 | 555.10 | 1353.04 | 1117.80 | 5.48 | 7.97 | 192.94 | 2222.71 | 8331.05 | 1.02 | 41.32 | 181.82 | 2387.87 | 8057.83 | 9.3030 | 0.02 | 333 | 2223 | 100.00 | 14.85 | 8.9326 | 125 |
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
# Initialize KMeans and StandardScaler
kmeans = KMeans(n_clusters=6, random_state=42)
scaler = StandardScaler()
# Apply clustering for specific datasets
if file_name in ['FD002.txt', 'FD004.txt']:
# Ensure required columns are present
if {'os1', 'os2', 'os3'}.issubset(df_train.columns):
# Select operating conditions and scale them
operating_conditions = df_train[['os1', 'os2', 'os3']]
scaled_conditions = scaler.fit_transform(operating_conditions)
# Apply k-means clustering and assign operation_mode
df_train['operation_mode'] = kmeans.fit_predict(scaled_conditions)
else:
raise ValueError("Columns 'os1', 'os2', and 'os3' are missing in the dataset!")
else:
print(f"No clustering applied for dataset {file_name}.")
if file_name in ['FD002.txt', 'FD004.txt']:
df_test['operation_mode'] = kmeans.predict(scaler.transform(df_test[['os1', 'os2', 'os3']]))
for df in [df_train, df_test]:
display(df.head())
| unit | time | os1 | os2 | os3 | s1 | s2 | s3 | s4 | s5 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s18 | s19 | s20 | s21 | rul | operation_mode | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 42.0049 | 0.8400 | 100.0 | 445.00 | 549.68 | 1343.43 | 1112.93 | 3.91 | 5.70 | 137.36 | 2211.86 | 8311.32 | 1.01 | 41.69 | 129.78 | 2387.99 | 8074.83 | 9.3335 | 0.02 | 330 | 2212 | 100.00 | 10.62 | 6.3670 | 125 | 5 |
| 1 | 1 | 2 | 20.0020 | 0.7002 | 100.0 | 491.19 | 606.07 | 1477.61 | 1237.50 | 9.35 | 13.61 | 332.10 | 2323.66 | 8713.60 | 1.07 | 43.94 | 312.59 | 2387.73 | 8046.13 | 9.1913 | 0.02 | 361 | 2324 | 100.00 | 24.37 | 14.6552 | 125 | 3 |
| 2 | 1 | 3 | 42.0038 | 0.8409 | 100.0 | 445.00 | 548.95 | 1343.12 | 1117.05 | 3.91 | 5.69 | 138.18 | 2211.92 | 8306.69 | 1.01 | 41.66 | 129.62 | 2387.97 | 8066.62 | 9.4007 | 0.02 | 329 | 2212 | 100.00 | 10.48 | 6.4213 | 125 | 5 |
| 3 | 1 | 4 | 42.0000 | 0.8400 | 100.0 | 445.00 | 548.70 | 1341.24 | 1118.03 | 3.91 | 5.70 | 137.98 | 2211.88 | 8312.35 | 1.02 | 41.68 | 129.80 | 2388.02 | 8076.05 | 9.3369 | 0.02 | 328 | 2212 | 100.00 | 10.54 | 6.4176 | 125 | 5 |
| 4 | 1 | 5 | 25.0063 | 0.6207 | 60.0 | 462.54 | 536.10 | 1255.23 | 1033.59 | 7.05 | 9.00 | 174.82 | 1915.22 | 7994.94 | 0.93 | 36.48 | 164.11 | 2028.08 | 7865.80 | 10.8366 | 0.02 | 305 | 1915 | 84.93 | 14.03 | 8.6754 | 125 | 2 |
| unit | time | os1 | os2 | os3 | s1 | s2 | s3 | s4 | s5 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s18 | s19 | s20 | s21 | rul | operation_mode | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 20.0072 | 0.7000 | 100.0 | 491.19 | 606.67 | 1481.04 | 1227.81 | 9.35 | 13.60 | 332.52 | 2323.67 | 8704.98 | 1.07 | 43.83 | 313.03 | 2387.78 | 8048.98 | 9.2229 | 0.02 | 362 | 2324 | 100.00 | 24.31 | 14.7007 | 125 | 3 |
| 1 | 1 | 2 | 24.9984 | 0.6200 | 60.0 | 462.54 | 536.22 | 1256.17 | 1031.48 | 7.05 | 9.00 | 174.46 | 1915.21 | 7999.94 | 0.93 | 36.11 | 163.61 | 2028.09 | 7863.46 | 10.8632 | 0.02 | 306 | 1915 | 84.93 | 14.36 | 8.5748 | 125 | 2 |
| 2 | 1 | 3 | 42.0000 | 0.8420 | 100.0 | 445.00 | 549.23 | 1340.13 | 1105.88 | 3.91 | 5.69 | 137.34 | 2211.93 | 8305.38 | 1.01 | 41.52 | 129.98 | 2387.95 | 8071.13 | 9.3960 | 0.02 | 328 | 2212 | 100.00 | 10.39 | 6.4365 | 125 | 5 |
| 3 | 1 | 4 | 42.0035 | 0.8402 | 100.0 | 445.00 | 549.19 | 1339.70 | 1107.26 | 3.91 | 5.69 | 137.23 | 2211.89 | 8301.00 | 1.01 | 41.73 | 129.48 | 2387.90 | 8078.89 | 9.3594 | 0.02 | 328 | 2212 | 100.00 | 10.56 | 6.2367 | 125 | 5 |
| 4 | 1 | 5 | 35.0079 | 0.8400 | 100.0 | 449.44 | 555.10 | 1353.04 | 1117.80 | 5.48 | 7.97 | 192.94 | 2222.71 | 8331.05 | 1.02 | 41.32 | 181.82 | 2387.87 | 8057.83 | 9.3030 | 0.02 | 333 | 2223 | 100.00 | 14.85 | 8.9326 | 125 | 0 |
Standard normalization (z-score normalization), the data is scaled to have a mean of 0 and a standard deviation of 1.
sensor_colums = ['s'+str(i) for i in range(1,22)]
print(sensor_colums)
['s1', 's2', 's3', 's4', 's5', 's6', 's7', 's8', 's9', 's10', 's11', 's12', 's13', 's14', 's15', 's16', 's17', 's18', 's19', 's20', 's21']
## Z-score Normalization
if file_name in ['FD001.txt', 'FD003.txt']:
# Dictionary to store mean and std for columns containing 's'
mean_std_dict = {}
# Calculate mean and std for 's' columns in df_train
for c in sensor_colums:
mean_std_dict[c + '_mean'] = df_train[c].mean()
mean_std_dict[c + '_std'] = df_train[c].std()
# Apply standard normalization to df_train
for c in sensor_colums:
df_train[c] = (df_train[c] - mean_std_dict[c + '_mean']) / (mean_std_dict[c + '_std'] + 1e-6)
# Apply standard normalization to df_test using df_train's mean and std
for c in sensor_colums:
df_test[c] = (df_test[c] - mean_std_dict[c + '_mean']) / (mean_std_dict[c + '_std'] + 1e-6)
if file_name in ['FD002.txt', 'FD004.txt']:
# Function to normalize each group based on its operation_mode
def z_score_normalize_by_mode(group):
group = group.copy() # To avoid modifying the original data
for c in sensor_colums:
group[c] = (group[c] - group[c].mean()) / (group[c].std() + 1e-6)
return group
# Normalize training data
df_train_normalized = df_train.groupby('operation_mode', group_keys=False).apply(z_score_normalize_by_mode)
# Normalize test data using train stats for each operation_mode
for mode in df_train['operation_mode'].unique():
train_group = df_train[df_train['operation_mode'] == mode]
test_group = df_test[df_test['operation_mode'] == mode]
for c in sensor_colums:
# Ensure the column dtype is compatible
df_test[c] = df_test[c].astype(float)
mean = train_group[c].mean()
std = train_group[c].std()
df_test.loc[df_test['operation_mode'] == mode, c] = (
df_test.loc[df_test['operation_mode'] == mode, c] - mean
) / (std+1e-6)
df_train = df_train_normalized.copy()
C:\Users\spiro\AppData\Local\Temp\ipykernel_9692\182794737.py:10: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.
df_train_normalized = df_train.groupby('operation_mode', group_keys=False).apply(z_score_normalize_by_mode)
if file_name in ['FD002.txt', 'FD004.txt']:
df_train = df_train.drop('operation_mode', axis = 1)
df_test = df_test.drop('operation_mode', axis = 1)
os_colums = ['os'+str(i) for i in range(1,4)]
print(os_colums)
mean_std_dict = {}
# Calculate mean and std for 's' columns in df_train
for c in os_colums:
mean_std_dict[c + '_mean'] = df_train[c].mean()
mean_std_dict[c + '_std'] = df_train[c].std()
# Apply standard normalization to df_train
for c in os_colums:
df_train[c] = (df_train[c] - mean_std_dict[c + '_mean']) / (mean_std_dict[c + '_std'] + 1e-6)
# Apply standard normalization to df_test using df_train's mean and std
for c in os_colums:
df_test[c] = (df_test[c] - mean_std_dict[c + '_mean']) / (mean_std_dict[c + '_std'] + 1e-6)
# Display the first few rows of both datasets
for df in [df_train, df_test]:
display(df.head())
['os1', 'os2', 'os3']
| unit | time | os1 | os2 | os3 | s1 | s2 | s3 | s4 | s5 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s18 | s19 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1.218323 | 0.864989 | 0.418618 | 0.000000e+00 | 0.432664 | -1.388852 | -1.286119 | 1.776357e-09 | -0.760397 | -1.517470 | -0.606409 | -0.878364 | -1.792634 | -1.211197 | -1.096189 | -0.384740 | -0.946674 | 0.020721 | 0.000000 | -0.330196 | 0.0 | 0.000000e+00 | -0.024111 | -0.079174 | 125 |
| 1 | 1 | 2 | -0.270016 | 0.415154 | 0.418618 | 1.136868e-07 | -2.477972 | -0.891902 | -1.361476 | 1.776357e-09 | -1.286994 | -1.471435 | -1.667478 | -0.880680 | -1.411996 | -1.422723 | -1.479263 | -1.914105 | -1.374161 | 0.009604 | -0.438767 | -2.282832 | 0.0 | 0.000000e+00 | -0.939500 | -0.633525 | 125 |
| 2 | 1 | 3 | 1.218249 | 0.867885 | 0.418618 | 0.000000e+00 | -1.139180 | -1.435904 | -0.828803 | 1.776357e-09 | -1.698290 | -0.724447 | -0.355415 | -1.128448 | -1.792634 | -1.315582 | -1.263988 | -0.462217 | -1.435947 | 0.952683 | 0.000000 | -0.933051 | 0.0 | 0.000000e+00 | -1.203853 | 0.682282 | 125 |
| 3 | 1 | 4 | 1.217992 | 0.864989 | 0.418618 | 0.000000e+00 | -1.677483 | -1.721248 | -0.720023 | 1.776357e-09 | -0.760397 | -0.917867 | -0.522745 | -0.822730 | -0.361048 | -1.245992 | -1.075214 | -0.268524 | -0.873968 | 0.067873 | 0.000000 | -1.535907 | 0.0 | 0.000000e+00 | -0.698249 | 0.630397 | 125 |
| 4 | 1 | 5 | 0.068489 | 0.159347 | -2.388771 | 5.684342e-08 | -1.658357 | -0.939885 | -1.768343 | 1.776357e-09 | -1.512372 | -0.911462 | -0.911903 | -1.444806 | -5.404283 | -0.866921 | -0.804348 | -1.032676 | -1.037340 | -0.383545 | 0.000000 | -1.278025 | 0.0 | -2.842171e-08 | -2.186010 | 1.153945 | 125 |
| unit | time | os1 | os2 | os3 | s1 | s2 | s3 | s4 | s5 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s18 | s19 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | -0.269665 | 0.414510 | 0.418618 | 1.136868e-07 | -1.296849 | -0.374676 | -2.451586 | 1.776357e-09 | -1.721582 | -1.273268 | -1.618666 | -1.364460 | -1.411996 | -1.811782 | -1.258376 | -1.676665 | -1.186464 | 0.495988 | -0.438767 | -1.687027 | 0.0 | 0.000000e+00 | -1.284064 | -0.196332 | 125 |
| 1 | 1 | 2 | 0.067955 | 0.157094 | -2.388771 | 5.684342e-08 | -1.347352 | -0.778231 | -2.050746 | 1.776357e-09 | -1.512372 | -1.257821 | -0.959334 | -1.068502 | -5.404283 | -2.402610 | -1.326402 | -0.987942 | -1.236464 | -0.018966 | 0.000000 | -0.579993 | 0.0 | -2.842171e-08 | 0.388681 | -0.162087 | 125 |
| 2 | 1 | 3 | 1.217992 | 0.871424 | 0.418618 | 0.000000e+00 | -0.536281 | -1.889723 | -2.068664 | 1.776357e-09 | -1.698290 | -1.536813 | -0.313583 | -1.199206 | -1.792634 | -1.802711 | -0.886439 | -0.539694 | -1.167174 | 0.887501 | 0.000000 | -1.535907 | 0.0 | 0.000000e+00 | -1.962258 | 0.895434 | 125 |
| 3 | 1 | 4 | 1.218228 | 0.865633 | 0.418618 | 0.000000e+00 | -0.622410 | -1.954988 | -1.915485 | 1.776357e-09 | -1.698290 | -1.643194 | -0.480912 | -1.435787 | -1.792634 | -1.072017 | -1.410813 | -0.733387 | -0.704719 | 0.379915 | 0.000000 | -1.535907 | 0.0 | 0.000000e+00 | -0.529715 | -1.906389 | 125 |
| 4 | 1 | 5 | 0.745026 | 0.864989 | 0.418618 | 1.705302e-07 | -0.959379 | -1.814420 | -1.210958 | -1.776357e-09 | -1.580330 | -1.437290 | -1.294164 | -1.399778 | -0.507686 | -2.020050 | -1.275919 | -0.976771 | -0.995511 | 0.208434 | 0.000000 | -0.450906 | 0.0 | 0.000000e+00 | -0.240601 | 0.035089 | 125 |
Below, we can see an example of all the sensor values of a specific engine in the training set (unit 4), as the engine progresses toward failure.
sample = 4
sample_df = df_train[df_train['unit'] == sample].copy()
# Select only sensor columns
n_sensors = len(sensor_colums)
rows = (n_sensors + 2) // 3 # Calculate rows for the grid
# Create subplots
fig, axes = plt.subplots(rows, 3, figsize=(15, 5 * rows))
# Flatten axes for easier iteration
axes = axes.flatten()
# Plot each sensor's data
for j, c in enumerate(sensor_colums):
axes[j].plot(sample_df['time'], sample_df[c], label='Sensor data')
axes[j].set_title(c)
# axes[j].set_xlabel('Time')
# axes[j].set_ylabel('Value')
axes[j].legend()
# Hide unused subplots
for ax in axes[n_sensors:]:
ax.axis('off')
# Adjust layout for better spacing
fig.tight_layout(pad=3.0)
plt.savefig(f'figures\\sensor_of_engine_{sample}_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
We plot the sensor data obtained from sensor 2 for 5 of the engines to show a general trend.
sensor = 's2'
for sample in range(1,6):
sample_df = df_train[df_train['unit'] == sample].copy()
sensordata = sample_df[sensor].to_numpy()
plt.plot(sensordata, label = "engine "+str(sample))
plt.grid()
plt.legend()
plt.ylabel('Sensor 2')
plt.xlabel('Cycle')
plt.savefig(f'figures\\sensor_{sensor}_of_engines_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
import seaborn as sns
import matplotlib.pyplot as plt
# Compute correlation matrix
corr_matrix = df_train.corr()
# Plot heatmap
plt.figure(figsize=(20, 14))
sns.heatmap(corr_matrix, annot=True, cmap="coolwarm", fmt=".2f")
plt.title("Feature Correlation Matrix")
plt.savefig(f'figures\\Feature_Correlation_Matrix_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
from sklearn.ensemble import RandomForestRegressor
import pandas as pd
Xt = df_train.iloc[:,2:-1]
yt = df_train.iloc[:,-1]
rf = RandomForestRegressor()
rf.fit(Xt, yt)
# Get feature importance
importance = pd.Series(rf.feature_importances_, index=Xt.columns)
importance = importance.sort_values(ascending=False)
print("Feature Importances:\n", importance)
Feature Importances: s11 0.587140 s9 0.095355 s14 0.082954 s4 0.031146 s15 0.027281 s6 0.025749 s13 0.021169 s8 0.015294 s3 0.015219 s12 0.015051 s2 0.015012 s7 0.013798 s21 0.012449 s20 0.011164 s17 0.009922 os1 0.009415 os2 0.006112 s10 0.003264 s16 0.001191 os3 0.000722 s1 0.000593 s5 0.000000 s18 0.000000 s19 0.000000 dtype: float64
# Combine thresholding and sorting
threshold = 1e-3
low_importance_features = importance[importance < threshold].sort_values()
print("Low importance features sorted:\n", low_importance_features)
Low importance features sorted: s5 0.000000 s18 0.000000 s19 0.000000 s1 0.000593 os3 0.000722 dtype: float64
We can notice here, that there are multiple sensors which are not changing its value. Perhaps, they are not useful features for prediction. Would they have similar behaviors for other engine units as well? Below show that yes (standard deviation is 0, or practically 0).
df_train[low_importance_features.index].describe()
| s5 | s18 | s19 | s1 | os3 | |
|---|---|---|---|---|---|
| count | 6.099400e+04 | 60994.0 | 6.099400e+04 | 6.099400e+04 | 6.099400e+04 |
| mean | 1.728769e-10 | 0.0 | -4.238047e-09 | 3.369189e-08 | -3.572867e-16 |
| std | 2.019544e-09 | 0.0 | 1.012389e-08 | 8.695509e-08 | 9.999999e-01 |
| min | -3.552714e-09 | 0.0 | -2.842171e-08 | -1.136868e-07 | -2.388771e+00 |
| 25% | -1.776357e-09 | 0.0 | 0.000000e+00 | 0.000000e+00 | 4.186184e-01 |
| 50% | 1.776357e-09 | 0.0 | 0.000000e+00 | 0.000000e+00 | 4.186184e-01 |
| 75% | 1.776357e-09 | 0.0 | 0.000000e+00 | 1.136868e-07 | 4.186184e-01 |
| max | 1.776357e-09 | 0.0 | 0.000000e+00 | 1.705302e-07 | 4.186184e-01 |
📌 Previously mentioned columns are dropped.
#Drop os3, s1, s5, s6, s10, s16, s18, s19 from both train and test
drop_cols1 = low_importance_features.index
df_train = df_train.drop(drop_cols1, axis = 1)
df_test = df_test.drop(drop_cols1, axis = 1)
for df in [df_train, df_test]:
display(df.head())
| unit | time | os1 | os2 | s2 | s3 | s4 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1.218323 | 0.864989 | 0.432664 | -1.388852 | -1.286119 | -0.760397 | -1.517470 | -0.606409 | -0.878364 | -1.792634 | -1.211197 | -1.096189 | -0.384740 | -0.946674 | 0.020721 | 0.000000 | -0.330196 | -0.024111 | -0.079174 | 125 |
| 1 | 1 | 2 | -0.270016 | 0.415154 | -2.477972 | -0.891902 | -1.361476 | -1.286994 | -1.471435 | -1.667478 | -0.880680 | -1.411996 | -1.422723 | -1.479263 | -1.914105 | -1.374161 | 0.009604 | -0.438767 | -2.282832 | -0.939500 | -0.633525 | 125 |
| 2 | 1 | 3 | 1.218249 | 0.867885 | -1.139180 | -1.435904 | -0.828803 | -1.698290 | -0.724447 | -0.355415 | -1.128448 | -1.792634 | -1.315582 | -1.263988 | -0.462217 | -1.435947 | 0.952683 | 0.000000 | -0.933051 | -1.203853 | 0.682282 | 125 |
| 3 | 1 | 4 | 1.217992 | 0.864989 | -1.677483 | -1.721248 | -0.720023 | -0.760397 | -0.917867 | -0.522745 | -0.822730 | -0.361048 | -1.245992 | -1.075214 | -0.268524 | -0.873968 | 0.067873 | 0.000000 | -1.535907 | -0.698249 | 0.630397 | 125 |
| 4 | 1 | 5 | 0.068489 | 0.159347 | -1.658357 | -0.939885 | -1.768343 | -1.512372 | -0.911462 | -0.911903 | -1.444806 | -5.404283 | -0.866921 | -0.804348 | -1.032676 | -1.037340 | -0.383545 | 0.000000 | -1.278025 | -2.186010 | 1.153945 | 125 |
| unit | time | os1 | os2 | s2 | s3 | s4 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | -0.269665 | 0.414510 | -1.296849 | -0.374676 | -2.451586 | -1.721582 | -1.273268 | -1.618666 | -1.364460 | -1.411996 | -1.811782 | -1.258376 | -1.676665 | -1.186464 | 0.495988 | -0.438767 | -1.687027 | -1.284064 | -0.196332 | 125 |
| 1 | 1 | 2 | 0.067955 | 0.157094 | -1.347352 | -0.778231 | -2.050746 | -1.512372 | -1.257821 | -0.959334 | -1.068502 | -5.404283 | -2.402610 | -1.326402 | -0.987942 | -1.236464 | -0.018966 | 0.000000 | -0.579993 | 0.388681 | -0.162087 | 125 |
| 2 | 1 | 3 | 1.217992 | 0.871424 | -0.536281 | -1.889723 | -2.068664 | -1.698290 | -1.536813 | -0.313583 | -1.199206 | -1.792634 | -1.802711 | -0.886439 | -0.539694 | -1.167174 | 0.887501 | 0.000000 | -1.535907 | -1.962258 | 0.895434 | 125 |
| 3 | 1 | 4 | 1.218228 | 0.865633 | -0.622410 | -1.954988 | -1.915485 | -1.698290 | -1.643194 | -0.480912 | -1.435787 | -1.792634 | -1.072017 | -1.410813 | -0.733387 | -0.704719 | 0.379915 | 0.000000 | -1.535907 | -0.529715 | -1.906389 | 125 |
| 4 | 1 | 5 | 0.745026 | 0.864989 | -0.959379 | -1.814420 | -1.210958 | -1.580330 | -1.437290 | -1.294164 | -1.399778 | -0.507686 | -2.020050 | -1.275919 | -0.976771 | -0.995511 | 0.208434 | 0.000000 | -0.450906 | -0.240601 | 0.035089 | 125 |
📌 Splitting Train and Validation Sets: Out of the 100 engines in the training set, I will randomly take out 20 engines for validation.
ratio = 0.8 #Ratio of training and validation datasets
units = np.arange(1, no_units+1)
no_selected = round(ratio * engine_numbers)
train_units = list(np.random.choice(units,no_selected, replace = False))
val_units = list(set(units) - set(train_units))
print(val_units)
train_data = df_train[df_train['unit'].isin(train_units)].copy()
val_data = df_train[df_train['unit'].isin(val_units)].copy()
[129, 5, 135, 136, 9, 139, 14, 15, 146, 24, 154, 28, 157, 33, 162, 163, 40, 41, 44, 45, 48, 53, 188, 190, 63, 191, 65, 71, 200, 202, 206, 208, 81, 82, 211, 213, 215, 88, 89, 217, 218, 221, 95, 228, 235, 237, 111, 242, 123, 124]
📌 The time series for sensor values were noisy. If the time-series values are $t_1, t_2, t_3, ..., t_n$, then the smoothed values $v_1, v_2, ..., v_n$ with the parameter $\beta$ follow the following formula:
$$v_0 = 0, v_{t} = \frac{\beta v_{t-1} + (1-\beta) x_{t}}{1-\beta^{t}}$$
#Smoothing Function: Exponentially Weighted Averages
def smooth(s, b = 0.98):
v = np.zeros(len(s)+1) #v_0 is already 0.
bc = np.zeros(len(s)+1)
for i in range(1, len(v)): #v_t = 0.95
v[i] = (b * v[i-1] + (1-b) * s[i-1])
bc[i] = 1 - b**i
sm = v[1:] / bc[1:]
return sm
# s = [1,2,3,4,5]
# print(s)
# print(f'After Smoothing: {smooth(s)}')
#Smoothing each time series for each engine in both training and test sets
beta = 0.98
# if file_name=='FD002' or file_name=='FD004':
# beta = 0.8
# Verify data integrity
assert 'unit' in df_train.columns, "The 'unit' column is missing in df_train"
assert 'unit' in df_test.columns, "The 'unit' column is missing in df_test"
def smooth_series(df, beta, unit_col='unit', sensor_prefix='s'):
"""
Smooths time series for each engine and sensor column in the dataset.
"""
for col in df.columns:
if sensor_prefix in col:
sm_list = []
for unit in df[unit_col].unique():
# Get sensor data for the current unit
unit_data = df[df[unit_col] == unit]
s = np.array(unit_data[col].copy())
sm = list(smooth(s, beta)) # Apply smoothing
sm_list.extend(sm) # Append smoothed data
# Check the length match
if len(sm_list) != len(df):
raise ValueError(f"Length mismatch for column {col}: sm_list={len(sm_list)}, df={len(df)}")
# Add smoothed column to DataFrame
df[col + '_smoothed'] = sm_list
return df
# Apply smoothing to training and test sets
df_train = smooth_series(df_train, beta)
df_test = smooth_series(df_test, beta)
Let's take a look at how smoothed values (salmon) look compared to the original series (lightblue), for a particular sensor (unit 10 from train)
df_train.columns
Index(['unit', 'time', 'os1', 'os2', 's2', 's3', 's4', 's6', 's7', 's8', 's9',
's10', 's11', 's12', 's13', 's14', 's15', 's16', 's17', 's20', 's21',
'rul', 'os1_smoothed', 'os2_smoothed', 's2_smoothed', 's3_smoothed',
's4_smoothed', 's6_smoothed', 's7_smoothed', 's8_smoothed',
's9_smoothed', 's10_smoothed', 's11_smoothed', 's12_smoothed',
's13_smoothed', 's14_smoothed', 's15_smoothed', 's16_smoothed',
's17_smoothed', 's20_smoothed', 's21_smoothed'],
dtype='object')
# Set scientific theme
sns.set_theme(style="whitegrid", font="Times New Roman")
# Select unit
sample_unit = 10
sample_df = df_train[df_train['unit'] == sample_unit].copy()
# Extract sensor columns
sensor_cols = [col for col in df_train.columns if col.startswith('s') and '_smoothed' not in col]
n_sensors = len(sensor_cols)
if n_sensors == 0:
raise ValueError("No valid sensor columns found in the dataset.")
# Define layout for compact visualization
num_cols = 5
num_rows = -(-n_sensors // num_cols)
fig, axes = plt.subplots(num_rows, num_cols, figsize=(9, 2.5 * num_rows), sharex=True, constrained_layout=True)
# Flatten axes
axes = axes.flatten()
# Define colors
colors = {"original": "tab:blue", "smoothed": "#D55E00", "line": "black"}
for j, col in enumerate(sensor_cols):
axes[j].plot(sample_df['time'], sample_df[col], color=colors["original"], label="Original", lw=0.8)
smoothed_col = f"{col}_smoothed"
if smoothed_col in sample_df.columns:
axes[j].plot(sample_df['time'], sample_df[smoothed_col], color=colors["smoothed"], label="Smoothed", lw=0.8)
# Grid and Labels
axes[j].set_title(col, fontsize=10)
axes[j].tick_params(axis='both', labelsize=8)
axes[j].grid(True, linestyle="--", linewidth=0.5, alpha=0.6)
# Hide unused subplots
for ax in axes[n_sensors:]:
ax.set_visible(False)
# Global legend for clarity
fig.legend(["Original", "Smoothed"], loc="upper center", ncol=2, bbox_to_anchor=(0.5, -0.02), fontsize=9, frameon=False)
# Save and display
output_dir = "figures"
os.makedirs(output_dir, exist_ok=True)
plt.savefig(os.path.join(output_dir, f"Smoothing_Sensors_{file_name}.png"), dpi=600, bbox_inches="tight", pad_inches=0.05)
plt.show()
#Remove the original series
for c in df_train.columns:
if ('s' in c) and ('smoothed' not in c):
df_train[c] = df_train[c+'_smoothed']
df_train.drop(c+'_smoothed', axis = 1, inplace = True)
for c in df_test.columns:
if ('s' in c) and ('smoothed' not in c):
df_test[c] = df_test[c+'_smoothed']
df_test.drop(c+'_smoothed', axis = 1, inplace = True)
for df in [df_train, df_test]:
display(df.head())
| unit | time | os1 | os2 | s2 | s3 | s4 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1.218323 | 0.864989 | 0.432664 | -1.388852 | -1.286119 | -0.760397 | -1.517470 | -0.606409 | -0.878364 | -1.792634 | -1.211197 | -1.096189 | -0.384740 | -0.946674 | 0.020721 | 0.000000 | -0.330196 | -0.024111 | -0.079174 | 125 |
| 1 | 1 | 2 | 0.466637 | 0.637800 | -1.037355 | -1.137867 | -1.324178 | -1.026355 | -1.494220 | -1.142303 | -0.879534 | -1.600393 | -1.318028 | -1.289661 | -1.157146 | -1.162577 | 0.015106 | -0.221599 | -1.316376 | -0.486428 | -0.359149 | 125 |
| 2 | 1 | 3 | 0.722252 | 0.716049 | -1.071985 | -1.239226 | -1.155706 | -1.254874 | -1.232428 | -0.874690 | -0.964187 | -1.665772 | -1.317196 | -1.280930 | -0.920808 | -1.255547 | 0.333967 | -0.146236 | -1.186011 | -0.730417 | -0.004969 | 125 |
| 3 | 1 | 4 | 0.849968 | 0.754420 | -1.227977 | -1.363408 | -1.043463 | -1.127483 | -1.151389 | -0.784020 | -0.927744 | -1.329641 | -1.298852 | -1.227932 | -0.752763 | -1.157242 | 0.265414 | -0.108562 | -1.276153 | -0.722130 | 0.158718 | 125 |
| 4 | 1 | 5 | 0.687294 | 0.630549 | -1.317566 | -1.275247 | -1.194355 | -1.207603 | -1.101445 | -0.810640 | -1.035376 | -2.177825 | -1.208941 | -1.139758 | -0.811030 | -1.132283 | 0.130326 | -0.085963 | -1.276543 | -1.026853 | 0.365886 | 125 |
| unit | time | os1 | os2 | s2 | s3 | s4 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | -0.269665 | 0.414510 | -1.296849 | -0.374676 | -2.451586 | -1.721582 | -1.273268 | -1.618666 | -1.364460 | -1.411996 | -1.811782 | -1.258376 | -1.676665 | -1.186464 | 0.495988 | -0.438767 | -1.687027 | -1.284064 | -0.196332 | 125 |
| 1 | 1 | 2 | -0.099150 | 0.284502 | -1.322355 | -0.578492 | -2.249142 | -1.615921 | -1.265467 | -1.285670 | -1.214986 | -3.428303 | -2.110180 | -1.292732 | -1.328825 | -1.211716 | 0.235910 | -0.217168 | -1.127919 | -0.439244 | -0.179037 | 125 |
| 2 | 1 | 3 | 0.348796 | 0.484108 | -1.055020 | -1.024428 | -2.187763 | -1.643934 | -1.357749 | -0.955073 | -1.209620 | -2.872028 | -2.005613 | -1.154556 | -1.060450 | -1.196568 | 0.457510 | -0.143311 | -1.266671 | -0.957205 | 0.186380 | 125 |
| 3 | 1 | 4 | 0.572785 | 0.582399 | -0.943568 | -1.264165 | -2.117617 | -1.657937 | -1.431287 | -0.832917 | -1.267886 | -2.593948 | -1.765094 | -1.220575 | -0.976190 | -1.069855 | 0.437519 | -0.106390 | -1.336033 | -0.847072 | -0.352772 | 125 |
| 4 | 1 | 5 | 0.608639 | 0.641224 | -0.946859 | -1.378707 | -1.928885 | -1.641782 | -1.432537 | -0.928931 | -1.295341 | -2.159668 | -1.818166 | -1.232095 | -0.976311 | -1.054380 | 0.389832 | -0.084244 | -1.151784 | -0.720828 | -0.272035 | 125 |
📌 When we look at the length of the trajectories for each unit in both the training and test sets, thus 200 engines, we see that for the training set that had the entire trajectories, have minimum length of 128. However, for the test set, we have a minimum value of 31. This means that if we were to predict the final RUL for every test engine unit, we cannot use a window size greater than 31 for training the model. Furthermore, if you take a look at the smoothed series visualization, I put a black vertical line for every graph, at about time step 10. This is because, when we smooth using exponentially weighted averages, seems like in the beginning, it seems a little shaky. Thus, I will be taking that part out of the training process. This will happen for test data as well. Thus, the maximum window size we can take is 31 - 10 = 21. I will go with 20.
print('training set time cycles:')
display(df_train.groupby('unit')['time'].max().describe())
print('test set time cycles:')
display(df_test.groupby('unit')['time'].max().describe())
training set time cycles:
count 248.000000 mean 245.943548 std 73.256391 min 128.000000 25% 189.250000 50% 233.000000 75% 291.250000 max 543.000000 Name: time, dtype: float64
test set time cycles:
count 248.000000 mean 166.185484 std 91.606522 min 19.000000 25% 102.000000 50% 153.500000 75% 218.500000 max 486.000000 Name: time, dtype: float64
n_features = len([c for c in df_train.columns if 's' in c]) + 1 #plus one for time
window = 15
print(f'number of features: {n_features}, window size: {window}')
number of features: 20, window size: 15
📌 Splitting Train and Validation Sets: Out of the 100 engines in the training set, I will randomly take out 20 engines for validation.
train_indices = list(train_data[(train_data['rul'] >= (window - 1)) & (train_data['time'] > 10)].index)
val_indices = list(val_data[(val_data['rul'] >= (window - 1)) & (val_data['time'] > 10)].index)
📌 Normalize the RUL index dividing by maximum value
rul_max = max(df_train['rul'])
df_train['rul'] = df_train['rul'] / rul_max
df_test['rul'] = df_test['rul'] / rul_max
rul_max
125
📌 Z-Normalize the time index
# max_time = max(df_train['time'])
df_train_mean = df_train['time'].mean()
df_train_std = df_train['time'].std()
df_train['time'] = (df_train['time'] - df_train_mean) / df_train_std
df_test['time'] = (df_test['time'] - df_train_mean) / df_train_std
display(df_train)
| unit | time | os1 | os2 | s2 | s3 | s4 | s6 | s7 | s8 | s9 | s10 | s11 | s12 | s13 | s14 | s15 | s16 | s17 | s20 | s21 | rul | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | -1.484104 | 1.218323 | 0.864989 | 0.432664 | -1.388852 | -1.286119 | -0.760397 | -1.517470 | -0.606409 | -0.878364 | -1.792634 | -1.211197 | -1.096189 | -0.384740 | -0.946674 | 0.020721 | 0.000000 | -0.330196 | -0.024111 | -0.079174 | 1.000 |
| 1 | 1 | -1.472974 | 0.466637 | 0.637800 | -1.037355 | -1.137867 | -1.324178 | -1.026355 | -1.494220 | -1.142303 | -0.879534 | -1.600393 | -1.318028 | -1.289661 | -1.157146 | -1.162577 | 0.015106 | -0.221599 | -1.316376 | -0.486428 | -0.359149 | 1.000 |
| 2 | 1 | -1.461843 | 0.722252 | 0.716049 | -1.071985 | -1.239226 | -1.155706 | -1.254874 | -1.232428 | -0.874690 | -0.964187 | -1.665772 | -1.317196 | -1.280930 | -0.920808 | -1.255547 | 0.333967 | -0.146236 | -1.186011 | -0.730417 | -0.004969 | 1.000 |
| 3 | 1 | -1.450713 | 0.849968 | 0.754420 | -1.227977 | -1.363408 | -1.043463 | -1.127483 | -1.151389 | -0.784020 | -0.927744 | -1.329641 | -1.298852 | -1.227932 | -0.752763 | -1.157242 | 0.265414 | -0.108562 | -1.276153 | -0.722130 | 0.158718 | 1.000 |
| 4 | 1 | -1.439583 | 0.687294 | 0.630549 | -1.317566 | -1.275247 | -1.194355 | -1.207603 | -1.101445 | -0.810640 | -1.035376 | -2.177825 | -1.208941 | -1.139758 | -0.811030 | -1.132283 | 0.130326 | -0.085963 | -1.276543 | -1.026853 | 0.365886 | 1.000 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 60989 | 248 | 0.508239 | 0.121028 | 0.152647 | 1.120924 | 0.971681 | 1.265227 | 0.799282 | -0.682815 | -0.501108 | -0.619780 | -0.269300 | 1.352772 | -0.621083 | -0.496962 | -0.932685 | 1.264394 | 0.276756 | 1.001991 | -0.945773 | -0.865997 | 0.032 |
| 60990 | 248 | 0.519369 | 0.085226 | 0.111781 | 1.155888 | 0.996001 | 1.290563 | 0.797572 | -0.691799 | -0.467899 | -0.626490 | -0.271036 | 1.374981 | -0.628321 | -0.461072 | -0.940386 | 1.280041 | 0.271074 | 1.020984 | -0.972959 | -0.886069 | 0.024 |
| 60991 | 248 | 0.530500 | 0.084878 | 0.112711 | 1.139622 | 1.001376 | 1.319797 | 0.800952 | -0.681105 | -0.485769 | -0.616503 | -0.262302 | 1.396260 | -0.637931 | -0.481980 | -0.943458 | 1.285953 | 0.265512 | 1.031103 | -0.964230 | -0.911317 | 0.016 |
| 60992 | 248 | 0.541630 | 0.108114 | 0.128231 | 1.167516 | 1.015638 | 1.340998 | 0.807401 | -0.684969 | -0.519986 | -0.616844 | -0.264327 | 1.408435 | -0.625176 | -0.520504 | -0.949622 | 1.296184 | 0.260067 | 1.040276 | -0.988154 | -0.921371 | 0.008 |
| 60993 | 248 | 0.552761 | 0.072631 | 0.087936 | 1.183535 | 1.036390 | 1.366365 | 0.805527 | -0.695433 | -0.482295 | -0.620102 | -0.266163 | 1.439875 | -0.630817 | -0.475905 | -0.959527 | 1.302955 | 0.254736 | 1.046569 | -1.007970 | -0.944476 | 0.000 |
60994 rows × 22 columns
📌 Prepare Training, Validation and Test Dataloaders.
For training process, I will take batches of 128.
For the validation process, I will take the entire validation set (all the windows avaiable).
For the test process, I will take the LAST window of each engine's given trajectory, thus it will have exactly 100 X(size 20 window)s and 100 y(RUL)s.
class data(Dataset):
def __init__(self, list_indices, df_train):
self.indices = list_indices
self.df_train = df_train
def __len__(self):
return len(self.indices)
def __getitem__(self, idx):
ind = self.indices[idx]
X_ = self.df_train.iloc[ind : ind + window, :].drop(['unit','rul'], axis = 1).copy().to_numpy()
y_ = self.df_train.iloc[ind + window - 1]['rul']
return X_, y_
train = data(train_indices, df_train)
val = data(val_indices, df_train)
trainloader = DataLoader(train, batch_size = 256, shuffle = True)
valloader = DataLoader(val, batch_size = 128, shuffle = True)
units = np.arange(1, no_units+1)
class test(Dataset):
def __init__(self, units, df_test):
self.units = units
self.df_test = df_test
def __len__(self):
return len(self.units)
def __getitem__(self, idx):
n = self.units[idx]
U = self.df_test[self.df_test['unit'] == n].copy()
X_ = U.reset_index().iloc[-window:,:].drop(['index', 'unit','rul'], axis = 1).copy().to_numpy()
y_ = U['rul'].min()
return X_, y_
test = test(units, df_test)
testloader = DataLoader(test, batch_size = 100)
dataiter = iter(trainloader)
x,y = next(dataiter)
x.shape
torch.Size([256, 15, 20])
2. Model Building and Training¶
📌 I will be using a Model incoprorating Attention Mechanism with pairwise interactions and positional awareness layers before the final output.
def anderson(f, x0, m=5, lam=1e-4, max_iter=50, tol=1e-2, beta=1.0, verbose=False):
"""Improved Anderson acceleration for fixed-point iteration."""
# Shape parameters
batch, channels, dim = x0.shape
# Storage for historical X and F
X = torch.zeros(batch, m, channels * dim, dtype=x0.dtype, device=x0.device)
F = torch.zeros(batch, m, channels * dim, dtype=x0.dtype, device=x0.device)
# Initialize X and F with the first two iterations
X[:, 0], F[:, 0] = x0.view(batch, -1), f(x0).view(batch, -1)
X[:, 1], F[:, 1] = F[:, 0], f(F[:, 0].view(batch, channels, dim)).view(batch, -1)
# Prepare the H matrix and y vector for Anderson acceleration
H = torch.zeros(batch, m + 1, m + 1, dtype=x0.dtype, device=x0.device)
H[:, 0, 1:] = H[:, 1:, 0] = 1 # First row and column for normalization
y = torch.zeros(batch, m + 1, 1, dtype=x0.dtype, device=x0.device)
y[:, 0, 0] = 1 # y vector for the Anderson step
res = [] # To store residuals
for k in range(2, max_iter):
n = min(k, m) # Number of iterations to consider
# Compute residual matrix G
G = F[:, :n] - X[:, :n]
# Compute the H matrix (Gramian + regularization)
GTG = torch.bmm(G, G.transpose(1, 2))
H[:, 1:n+1, 1:n+1] = GTG + lam * torch.eye(n, dtype=x0.dtype, device=x0.device)[None]
# Solve for alpha using least squares
try:
alpha = torch.linalg.solve(H[:, :n+1, :n+1], y[:, :n+1])[:, 1:n+1, 0] # (batch x n)
except RuntimeError as e:
if verbose:
print(f"Solver failed at iteration {k}: {e}")
break
# Update X and F
X_update = beta * (alpha[:, None] @ F[:, :n])[:, 0] + (1 - beta) * (alpha[:, None] @ X[:, :n])[:, 0]
idx = k % m
X[:, idx] = X_update
F[:, idx] = f(X[:, idx].view(batch, channels, dim)).view(batch, -1)
# Compute residuals
residual = (F[:, idx] - X[:, idx]).norm(dim=1) / (1e-5 + X[:, idx].norm(dim=1).clamp(min=1e-5))
res.append(residual.mean().item()) # Mean residual across batch
if verbose:
print(f"Iteration {k}, residual: {res[-1]:.6f}")
# Check for convergence
if res[-1] < tol:
if verbose:
print(f"Converged at iteration {k}, residual: {res[-1]:.6f}")
break
return X[:, idx].view(batch, channels, dim), (res, k)
import torch.autograd as autograd
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print('Using device:', device)
class DEQFixedPoint(nn.Module):
def __init__(self, f, solver, **kwargs):
super().__init__()
self.f = f
self.solver = solver
self.kwargs = kwargs
def forward(self, x):
# Compute forward pass and re-engage autograd tape
b, c, f = x.shape
with torch.no_grad(): # Ensure no gradients for intermediate calculations
z, self.stats = self.solver(lambda z: self.f(z, x), torch.zeros((b, 4 * c, f), device=device), **self.kwargs)
z = self.f(z, x)
# Set z to require gradients before the hook
z = z.requires_grad_() # Ensure z requires gradients
# Set up Jacobian vector product (without additional forward calls)
z0 = z.clone().detach().requires_grad_() # Ensure z0 requires gradients
f0 = self.f(z0, x)
def backward_hook(grad):
# Perform the backward pass using the solver
g, self.backward_res = self.solver(
lambda y: autograd.grad(f0, z0, y, retain_graph=True)[0] + grad,
grad, **self.kwargs)
return g
# Register the backward hook
z.register_hook(backward_hook)
return z, self.stats
Using device: cuda
import torch.optim.lr_scheduler as lr_scheduler
# class ScaledLeakyReLU(nn.Module):
# def __init__(self, negative_slope=0.1, scaled_slope=0.9):
# super(ScaledLeakyReLU, self).__init__()
# self.negative_slope = negative_slope
# self.scaled_slope = scaled_slope
# def forward(self, x):
# return torch.where(x > 0, self.scaled_slope * x, self.negative_slope * x)
class DualInputAttention(nn.Module):
def __init__(self, channels, feature_size):
super(DualInputAttention, self).__init__()
self.feature_size = feature_size
self.query_proj = nn.Linear(self.feature_size, self.feature_size)
self.key_proj = nn.Linear(self.feature_size, self.feature_size)
self.value_proj = nn.Linear(self.feature_size, self.feature_size)
self.output_proj = nn.Linear(self.feature_size, self.feature_size)
self.query_proj.weight.data.normal_(0, 0.01)
self.key_proj.weight.data.normal_(0, 0.01)
self.value_proj.weight.data.normal_(0, 0.01)
self.output_proj.weight.data.normal_(0, 0.01)
def forward(self, input1, input2):
"""
Args:
input1: Tensor of shape (batch, channels, feature_size)
input2: Tensor of shape (batch, channels, feature_size)
Returns:
Tensor of shape (batch, channels, feature_size), attended features.
"""
# Project inputs
query = self.query_proj(input1) # (batch, channels, feature_size)
key = self.key_proj(input2) # (batch, channels, feature_size)
value = self.value_proj(input2) # (batch, channels, feature_size)
# Compute similarity (scaled dot product)
scores = torch.einsum('bci,bcj->bcij', query, key) # (batch, channels, feature_size, feature_size)
scores = scores / (self.feature_size ** 0.5)
# Compute attention weights
attention_weights = F.softmax(scores, dim=-1) # (batch, channels, feature_size, feature_size)
# Compute attended values
attended = torch.einsum('bcij,bcj->bci', attention_weights, value) # (batch, channels, feature_size)
# Combine with input1
output = self.output_proj(attended + input1) # (batch, channels, feature_size)
return output
class DELayer(nn.Module):
def __init__(self, n_channels, feature_size, kernel_size=3):
super().__init__()
num_groups = 4
self.conv0 = nn.Conv1d(n_channels, 2 * n_channels, kernel_size=1, padding='same', bias=False)
self.conv1 = nn.Conv1d(2 * n_channels, 2 * n_channels, kernel_size=kernel_size, padding='same', bias=False)
self.conv2 = nn.Conv1d(2 * n_channels, 2 * n_channels, kernel_size=kernel_size, padding='same', bias=False)
self.norm1 = nn.GroupNorm(num_groups, 4 * n_channels)
self.norm2 = nn.GroupNorm(num_groups, 4 * n_channels)
# self.norm3 = nn.GroupNorm(num_groups, 2 * n_channels)
# self.scaled_leaky_relu = ScaledLeakyReLU(negative_slope=0.1, scaled_slope=0.9)
self.attention = DualInputAttention(4 * n_channels, feature_size)
self.conv0.weight.data.normal_(0, 0.01)
self.conv1.weight.data.normal_(0, 0.01)
self.conv2.weight.data.normal_(0, 0.01)
def forward(self, z0, x):
# z = self.conv2(self.norm1(F.relu(self.conv1(z0))))
# out = self.attention (self.norm2(z),self.norm3(x))
# z = self.norm4(F.relu(out))
x = F.relu(self.conv0(x))
x_conv1 = F.relu(self.conv1(x))
x_conv2 = F.relu(self.conv2(x_conv1))
x_concat = torch.cat([x_conv1, x_conv2], dim=1)
x = self.norm1(x_concat)
# z = self.norm2(z0)
z = self.norm2(F.relu(self.attention(z0, x)))
return z
# Define the Bayesian Neural Network with Dropout
class BNN(nn.Module):
def __init__(self, initial_channels, hidden_dim, dropout):
super(BNN, self).__init__()
self.dropout_value = dropout
self.fc1 = nn.Linear(4 * initial_channels, 8* hidden_dim)
self.dropout1 = nn.Dropout(self.dropout_value) # Dropout for Monte Carlo Dropout
self.fc2 = nn.Linear(8 * hidden_dim, 4 * hidden_dim)
self.dropout2 = nn.Dropout(self.dropout_value)
self.fc3 = nn.Linear(4 * hidden_dim, 1)
def forward(self, x, training=True):
x = torch.relu(self.fc1(x))
x = self.dropout1(x) if training else x # Apply dropout during training and inference
x = torch.relu(self.fc2(x))
x = self.dropout2(x) if training else x # Apply dropout during training and inference
return self.fc3(x)
class RegressorModel(nn.Module):
def __init__(self, initial_channels, feature_size, hidden_dim, dropout):
super(RegressorModel, self).__init__()
# Layers
# self.conv =nn.Conv1d(initial_channels, 2 * initial_channels, kernel_size=1, bias=False)
# self.norm = nn.GroupNorm(4, 2 * initial_channels)
self.f = DELayer(n_channels= initial_channels, feature_size=feature_size )
self.DEQ = DEQFixedPoint(self.f, anderson, tol=1e-4, max_iter=200, beta=1.0)
self.pool = nn.AdaptiveAvgPool1d(1)
self.feedforward = BNN(initial_channels, hidden_dim, dropout)
def forward(self, x, training=True):
# x = self.norm(F.relu(self.conv(x)))
x = self.DEQ(x)
x = self.pool(x[0]).squeeze(-1) # Apply Deep Equilibrium Model
# x = torch.flatten(x, start_dim=1)
output = self.feedforward(x, training) # Apply feedforward Model for regression
return output
learning_rate = 1e-3
init_channels = window
inner_channels = window
n_hidden_units = window
set_dropout = 0.4
model = RegressorModel(initial_channels=init_channels, feature_size=n_features, hidden_dim=n_hidden_units, dropout=set_dropout ).to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)
scheduler = lr_scheduler.StepLR(optimizer, step_size=10, gamma=0.5)
X= torch.randn((64, window, n_features), device=device)
f = DELayer(n_channels=init_channels, feature_size= n_features).to(device)
DEQ = DEQFixedPoint(f, anderson, tol=1e-4, max_iter=100, beta=0.9).to(device)
out = DEQ(X)[0]
print(out.shape)
out2 = nn.AdaptiveAvgPool1d(1)(out).squeeze(-1)
print(out2.shape)
out3 = model(X)[0]
print(out3.shape)
torch.Size([64, 60, 20]) torch.Size([64, 60]) torch.Size([1])
📌 Validation and Test Functions
# def validation(loss_fn, num_samples=100):
# model.train() # Enable stochastic behavior (dropout during inference)
# total_loss = 0.0
# total_samples = 0
# with torch.no_grad():
# for X, y in valloader:
# X, y = X.to(device).float(), y.to(device).float()
# stochastic_predictions = []
# # Perform stochastic forward passes
# for _ in range(num_samples):
# y_pred = model(X, training=True).squeeze()
# stochastic_predictions.append(y_pred.cpu().numpy())
# # Convert to NumPy array
# stochastic_predictions = np.array(stochastic_predictions) # Shape: (num_samples, batch_size)
# # Calculate mean prediction
# mean_pred = torch.tensor(np.mean(stochastic_predictions, axis=0)).to(device)
# # Compute batch loss
# loss = loss_fn(mean_pred, y)
# total_loss += loss.item() * X.size(0) # Accumulate the weighted batch loss
# total_samples += X.size(0)
# avg_val_loss = total_loss / total_samples # Calculate average loss
# return avg_val_loss
def validation(loss_fn):
model.eval() # Set the model to evaluation mode
total_loss = 0.0
total_samples = 0
with torch.no_grad(): # Disable gradient computation for validation
for X, y in valloader: # Iterate over all batches in the validation loader
X, y = X.to(device).to(torch.float32), y.to(device).to(torch.float32)
# Use enable_grad for specific models if required
if isinstance(model, RegressorModel): # Replace with the DEQ model class
with torch.enable_grad():
y_pred = model(X).squeeze()
else:
y_pred = model(X).squeeze()
y = y.squeeze()
loss = loss_fn(y_pred, y)
total_loss += loss.item() * X.size(0) # Accumulate the weighted batch loss
total_samples += X.size(0) # Count the total samples processed
avg_val_loss = total_loss / total_samples # Calculate average loss
return avg_val_loss
loss_L1 = nn.L1Loss()
def test():
model.train() # Enable stochastic behavior (dropout during inference)
total_loss_MSE = 0.0
total_loss_L1 = 0.0
total_ASUE = 0.0
all_stochastic_predictions = []
all_pred_values = []
all_true_values = []
all_uncertainties = [] # Collect uncertainties
num_samples = 100
with torch.no_grad():
for X, y in testloader:
X, y = X.to(device).float(), y.to(device).float()
stochastic_predictions = []
# Perform stochastic forward passes
for _ in range(num_samples):
y_pred = model(X, training=True).squeeze()
stochastic_predictions.append(y_pred.cpu().numpy())
# Convert to NumPy array
stochastic_predictions = np.array(stochastic_predictions) # Shape: (num_samples, batch_size)
# Calculate mean and uncertainty
mean_pred = torch.tensor(np.mean(stochastic_predictions, axis=0)).to(device)
uncertainty_batch = np.var(stochastic_predictions, axis=0) # Variance as uncertainty
# Compute batch losses
loss_MSE = torch.mean((mean_pred - y) ** 2).item()
loss_L1_val = loss_L1(mean_pred, y).item()
ASUE = torch.mean(torch.relu(y - mean_pred)).item()
# Aggregate results
total_loss_MSE += loss_MSE * X.size(0)
total_loss_L1 += loss_L1_val * X.size(0)
total_ASUE += ASUE * X.size(0)
# Collect predictions, true values, and uncertainties
all_stochastic_predictions.append(stochastic_predictions)
all_pred_values.append(mean_pred.cpu().numpy())
all_true_values.append(y.cpu().numpy())
all_uncertainties.append(uncertainty_batch)
# Average metrics
total_samples = len(testloader.dataset)
avg_loss_MSE = total_loss_MSE / total_samples
avg_loss_L1 = total_loss_L1 / total_samples
avg_ASUE = total_ASUE / total_samples
# Concatenate all predictions, uncertainties, and true values
all_pred_values = np.concatenate([x.flatten() for x in all_pred_values], axis=0)
all_stochastic_predictions = np.concatenate(all_stochastic_predictions, axis=1)
all_true_values = np.concatenate(all_true_values, axis=0)
all_uncertainties = np.concatenate(all_uncertainties, axis=0) # Ensure consistency
return avg_loss_MSE, avg_loss_L1, avg_ASUE, all_uncertainties, all_stochastic_predictions, all_pred_values, all_true_values
📌 Training Loop: I have trained using Adam Optimizer for 35 epochs with learning rate = 0.01 and scheduler with degration 0.5 after 8 epochs
# Initialize lists for tracking losses
T, V = [], []
epochs = 35
loss_fn = nn.MSELoss()
siterations, sresiduals = [], []
stats = {}
# Define forward hook for DEQ statistics
def forward_hook(module, input, output):
stats[module] = output
# Register forward hook
hook_handle = model.DEQ.register_forward_hook(forward_hook)
best_val_loss = float('inf') # Track best validation loss
pbar = tqdm(range(epochs), desc="Training Progress", dynamic_ncols=True)
for epoch in pbar:
model.train()
epoch_loss = 0
sresiduals_list, siterations_list = [], []
for X, y in trainloader:
X, y = X.to(device, dtype=torch.float32), y.to(device, dtype=torch.float32)
# Forward pass
y_pred = model(X).squeeze()
y = y.squeeze()
loss = loss_fn(y_pred, y)
epoch_loss += loss.item()
# Backward pass
optimizer.zero_grad()
loss.backward()
optimizer.step()
# Collect DEQ statistics
iterations = stats.get(model.DEQ, (None, [0, 0]))[1]
if iterations and len(iterations) >= 2:
sresiduals_list.append(iterations[0][-1])
siterations_list.append(iterations[1])
# Step the scheduler
scheduler.step()
# Validation loss
val_loss = validation(loss_fn)
T.append(epoch_loss / len(trainloader))
V.append(val_loss)
# Update tqdm bar
pbar.set_postfix({'Train Loss': T[-1], 'Val Loss': val_loss})
# Track DEQ statistics
sresiduals.append(np.mean(sresiduals_list) if sresiduals_list else 0)
siterations.append(np.mean(siterations_list) if siterations_list else 0)
# Save best model
if val_loss < best_val_loss:
best_val_loss = val_loss
torch.save(model.state_dict(), 'best_model.pth')
print(f'Epoch {epoch + 1}/{epochs} - Train Loss: {T[-1]:.4f}, Val Loss: {V[-1]:.4f}')
# Clean up hook handle
hook_handle.remove()
Training Progress: 3%|▎ | 1/35 [01:23<47:19, 83.52s/it, Train Loss=0.0979, Val Loss=0.0243]
Epoch 1/35 - Train Loss: 0.0979, Val Loss: 0.0243
Training Progress: 6%|▌ | 2/35 [02:48<46:26, 84.44s/it, Train Loss=0.033, Val Loss=0.0152]
Epoch 2/35 - Train Loss: 0.0330, Val Loss: 0.0152
Training Progress: 9%|▊ | 3/35 [04:07<43:48, 82.13s/it, Train Loss=0.0237, Val Loss=0.0131]
Epoch 3/35 - Train Loss: 0.0237, Val Loss: 0.0131
Training Progress: 11%|█▏ | 4/35 [05:31<42:40, 82.60s/it, Train Loss=0.021, Val Loss=0.0121]
Epoch 4/35 - Train Loss: 0.0210, Val Loss: 0.0121
Training Progress: 14%|█▍ | 5/35 [06:55<41:38, 83.27s/it, Train Loss=0.0192, Val Loss=0.0121]
Epoch 5/35 - Train Loss: 0.0192, Val Loss: 0.0121
Training Progress: 17%|█▋ | 6/35 [08:26<41:32, 85.94s/it, Train Loss=0.0175, Val Loss=0.0124]
Epoch 6/35 - Train Loss: 0.0175, Val Loss: 0.0124
Training Progress: 20%|██ | 7/35 [09:57<40:48, 87.43s/it, Train Loss=0.0164, Val Loss=0.0117]
Epoch 7/35 - Train Loss: 0.0164, Val Loss: 0.0117
Training Progress: 23%|██▎ | 8/35 [11:26<39:32, 87.88s/it, Train Loss=0.0152, Val Loss=0.0118]
Epoch 8/35 - Train Loss: 0.0152, Val Loss: 0.0118
Training Progress: 26%|██▌ | 9/35 [12:51<37:46, 87.16s/it, Train Loss=0.0145, Val Loss=0.0118]
Epoch 9/35 - Train Loss: 0.0145, Val Loss: 0.0118
Training Progress: 29%|██▊ | 10/35 [14:17<36:07, 86.69s/it, Train Loss=0.0137, Val Loss=0.0112]
Epoch 10/35 - Train Loss: 0.0137, Val Loss: 0.0112
Training Progress: 31%|███▏ | 11/35 [15:41<34:21, 85.89s/it, Train Loss=0.0127, Val Loss=0.0111]
Epoch 11/35 - Train Loss: 0.0127, Val Loss: 0.0111
Training Progress: 34%|███▍ | 12/35 [17:06<32:49, 85.62s/it, Train Loss=0.0124, Val Loss=0.011]
Epoch 12/35 - Train Loss: 0.0124, Val Loss: 0.0110
Training Progress: 37%|███▋ | 13/35 [18:30<31:13, 85.16s/it, Train Loss=0.0121, Val Loss=0.0114]
Epoch 13/35 - Train Loss: 0.0121, Val Loss: 0.0114
Training Progress: 40%|████ | 14/35 [19:53<29:34, 84.48s/it, Train Loss=0.0117, Val Loss=0.0117]
Epoch 14/35 - Train Loss: 0.0117, Val Loss: 0.0117
Training Progress: 43%|████▎ | 15/35 [21:16<28:02, 84.13s/it, Train Loss=0.0116, Val Loss=0.0114]
Epoch 15/35 - Train Loss: 0.0116, Val Loss: 0.0114
Training Progress: 46%|████▌ | 16/35 [22:41<26:42, 84.36s/it, Train Loss=0.0111, Val Loss=0.0114]
Epoch 16/35 - Train Loss: 0.0111, Val Loss: 0.0114
Training Progress: 49%|████▊ | 17/35 [24:10<25:40, 85.59s/it, Train Loss=0.011, Val Loss=0.0114]
Epoch 17/35 - Train Loss: 0.0110, Val Loss: 0.0114
Training Progress: 51%|█████▏ | 18/35 [25:36<24:16, 85.69s/it, Train Loss=0.0107, Val Loss=0.0114]
Epoch 18/35 - Train Loss: 0.0107, Val Loss: 0.0114
Training Progress: 54%|█████▍ | 19/35 [27:03<23:01, 86.31s/it, Train Loss=0.0104, Val Loss=0.0129]
Epoch 19/35 - Train Loss: 0.0104, Val Loss: 0.0129
Training Progress: 57%|█████▋ | 20/35 [28:30<21:36, 86.45s/it, Train Loss=0.0102, Val Loss=0.0118]
Epoch 20/35 - Train Loss: 0.0102, Val Loss: 0.0118
Training Progress: 60%|██████ | 21/35 [29:57<20:10, 86.50s/it, Train Loss=0.00975, Val Loss=0.0112]
Epoch 21/35 - Train Loss: 0.0097, Val Loss: 0.0112
Training Progress: 63%|██████▎ | 22/35 [31:23<18:43, 86.39s/it, Train Loss=0.0095, Val Loss=0.0114]
Epoch 22/35 - Train Loss: 0.0095, Val Loss: 0.0114
Training Progress: 66%|██████▌ | 23/35 [32:53<17:30, 87.54s/it, Train Loss=0.00955, Val Loss=0.0113]
Epoch 23/35 - Train Loss: 0.0096, Val Loss: 0.0113
Training Progress: 69%|██████▊ | 24/35 [34:22<16:05, 87.81s/it, Train Loss=0.0094, Val Loss=0.0114]
Epoch 24/35 - Train Loss: 0.0094, Val Loss: 0.0114
Training Progress: 71%|███████▏ | 25/35 [35:49<14:36, 87.64s/it, Train Loss=0.00927, Val Loss=0.0122]
Epoch 25/35 - Train Loss: 0.0093, Val Loss: 0.0122
Training Progress: 74%|███████▍ | 26/35 [37:19<13:14, 88.29s/it, Train Loss=0.00907, Val Loss=0.0119]
Epoch 26/35 - Train Loss: 0.0091, Val Loss: 0.0119
Training Progress: 77%|███████▋ | 27/35 [38:47<11:45, 88.19s/it, Train Loss=0.00903, Val Loss=0.0128]
Epoch 27/35 - Train Loss: 0.0090, Val Loss: 0.0128
Training Progress: 80%|████████ | 28/35 [40:15<10:18, 88.41s/it, Train Loss=0.00891, Val Loss=0.0129]
Epoch 28/35 - Train Loss: 0.0089, Val Loss: 0.0129
Training Progress: 83%|████████▎ | 29/35 [41:45<08:52, 88.71s/it, Train Loss=0.00868, Val Loss=0.0114]
Epoch 29/35 - Train Loss: 0.0087, Val Loss: 0.0114
Training Progress: 86%|████████▌ | 30/35 [43:19<07:31, 90.26s/it, Train Loss=0.00862, Val Loss=0.0117]
Epoch 30/35 - Train Loss: 0.0086, Val Loss: 0.0117
Training Progress: 89%|████████▊ | 31/35 [44:53<06:05, 91.48s/it, Train Loss=0.00843, Val Loss=0.0115]
Epoch 31/35 - Train Loss: 0.0084, Val Loss: 0.0115
Training Progress: 91%|█████████▏| 32/35 [46:26<04:35, 91.89s/it, Train Loss=0.00831, Val Loss=0.012]
Epoch 32/35 - Train Loss: 0.0083, Val Loss: 0.0120
Training Progress: 94%|█████████▍| 33/35 [48:00<03:05, 92.58s/it, Train Loss=0.00828, Val Loss=0.0118]
Epoch 33/35 - Train Loss: 0.0083, Val Loss: 0.0118
Training Progress: 97%|█████████▋| 34/35 [49:33<01:32, 92.54s/it, Train Loss=0.00816, Val Loss=0.0124]
Epoch 34/35 - Train Loss: 0.0082, Val Loss: 0.0124
Training Progress: 100%|██████████| 35/35 [51:08<00:00, 87.66s/it, Train Loss=0.00827, Val Loss=0.0117]
Epoch 35/35 - Train Loss: 0.0083, Val Loss: 0.0117
print(sresiduals)
print(siterations)
[4.629765501090862e-05, 6.0786133999076704e-05, 5.565338161368557e-05, 4.768837594425815e-05, 6.850324595244845e-05, 6.849385082987706e-05, 7.143836879803594e-05, 7.838040940885887e-05, 6.501814759969927e-05, 6.24514488113481e-05, 6.398383843483062e-05, 7.08309207131902e-05, 7.270502892213192e-05, 6.984984097365415e-05, 7.391439238174001e-05, 8.011613920453508e-05, 8.114433523249233e-05, 8.302382069137652e-05, 8.456120009890751e-05, 7.99294200480163e-05, 7.098811853888206e-05, 7.852117562643523e-05, 8.413747449129344e-05, 8.404635754718477e-05, 8.78115205163932e-05, 7.959117566886477e-05, 8.363073995793075e-05, 8.344922634721517e-05, 8.170895713657617e-05, 8.419091284659496e-05, 8.902548151402084e-05, 9.045941484611492e-05, 8.287076220873386e-05, 8.499276242049703e-05, 7.862451677156483e-05] [9.029069767441861, 9.034883720930232, 8.063953488372093, 8.0, 8.36046511627907, 9.372093023255815, 10.412790697674419, 10.505813953488373, 9.88953488372093, 9.575581395348838, 9.69186046511628, 9.633720930232558, 9.186046511627907, 9.180232558139535, 9.459302325581396, 9.44186046511628, 9.738372093023257, 10.337209302325581, 11.029069767441861, 11.470930232558139, 11.953488372093023, 12.0, 12.05813953488372, 12.05232558139535, 12.13953488372093, 12.093023255813954, 12.203488372093023, 12.372093023255815, 12.75, 12.587209302325581, 13.075581395348838, 13.05232558139535, 13.540697674418604, 13.261627906976743, 13.784883720930232]
torch.save(model.state_dict(), f'saved_models\\DEM_{file_name}')
# Perform testing with Monte Carlo Dropout
mse, l1, asue, uncertainty, all_y_pred, y_pred, y = test()
# Calculate RMSE and scale by RUL max value
rmse = rul_max * np.sqrt(mse)
# Print results
print(f"Test Results:\n"
f"RMSE: {round(rmse, 2)}\n"
f"L1 Loss: {round(l1, 2)}\n"
f"ASUE: {round(asue, 2)}")
# Mean uncertainty
mean_uncertainty = np.mean(uncertainty)
print(f"Mean Uncertainty: {round(mean_uncertainty, 4)}")
# Optional: Visualize predictions, true values, and uncertainties
import matplotlib.pyplot as plt
plt.figure(figsize=(14, 8))
# True vs predicted RUL
plt.plot(range(len(y)), y, label="True Values", c="lightseagreen", marker=".", alpha=0.7)
plt.plot(range(len(y_pred)), y_pred, label="Predicted Values", c="salmon", marker=".", alpha=0.7)
# Uncertainty as error bars
plt.fill_between(
range(len(y_pred)),
y_pred - np.sqrt(uncertainty),
y_pred + np.sqrt(uncertainty),
color="salmon", alpha=0.2, label="Uncertainty"
)
# Add labels and legend
plt.xlabel("Sample Index", fontsize=14)
plt.ylabel("RUL", fontsize=14)
plt.title("True vs Predicted RUL with Uncertainty", fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.savefig(f'figures\\Responses_with_Uncertainty_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
Test Results: RMSE: 11.79 L1 Loss: 0.07 ASUE: 0.03 Mean Uncertainty: 0.0034000000450760126
print(y)
print(y_pred)
print(1*(y_pred>y))
print(all_y_pred.shape)
[0.176 0.312 0.856 0.6 1. 0.624 0.752 0.112 0.792 1. 1. 0.056 0.568 0.84 0.096 1. 1. 0.832 1. 0.656 0.728 0.088 0.208 1. 0.312 0.736 0.608 0.992 0.512 0.944 0.048 0.176 1. 1. 0.288 0.584 0.712 0.088 1. 0.08 0.776 0.24 0.336 0.48 0.68 1. 0.272 0.36 0.192 0.688 0.952 1. 1. 1. 1. 0.536 0.776 0.064 1. 1. 0.408 0.264 1. 0.368 0.096 1. 0.368 0.368 0.096 0.264 0.12 1. 0.184 0.712 0.992 1. 0.2 0.592 0.624 0.912 0.768 0.08 1. 1. 0.92 0.56 0.752 0.448 0.688 0.768 0.4 0.584 1. 1. 1. 0.568 0.84 0.904 0.296 0.056 0.104 0.176 0.072 0.96 0.8 0.856 0.328 1. 1. 0.472 0.144 0.528 0.104 0.112 1. 0.104 0.6 0.064 0.872 1. 0.328 1. 0.184 0.688 1. 0.12 1. 1. 0.96 1. 0.808 0.928 1. 0.288 0.056 0.976 1. 0.704 1. 1. 1. 0.864 0.424 1. 0.472 0.8 0.448 1. 0.608 0.456 0.248 0.704 1. 0.272 0.056 1. 1. 0.048 0.176 0.664 0.656 0.672 0.76 1. 0.888 0.576 0.872 0.696 1. 1. 1. 0.096 0.064 0.08 0.984 0.824 0.096 0.848 0.096 0.256 0.296 0.928 0.12 0.08 0.368 1. 0.192 1. 0.448 0.344 1. 0.568 0.832 0.12 1. 0.712 0.288 0.088 0.736 0.768 0.472 0.104 1. 1. 1. 0.872 0.928 0.728 0.088 0.704 0.864 0.608 0.112 0.712 1. 0.136 0.528 1. 0.328 1. 0.584 0.312 0.464 0.112 1. 0.704 1. 1. 0.96 0.784 0.264 1. 0.88 0.544 0.192 0.6 0.144 0.128 1. 0.784 1. 0.648 0.944 0.28 1. 1. 0.896 0.208] [0.2462424 0.22612968 0.9726383 0.74565536 0.98987794 0.5629638 0.9876078 0.14675337 0.77300537 0.97599757 0.97932214 0.06907845 0.60648227 0.92693496 0.13410877 1.0029427 0.94054323 0.74515706 0.98858225 0.68648803 0.6579274 0.09368446 0.22357674 0.9828942 0.42533615 0.9052162 0.5663383 0.9802998 0.5275922 0.9079585 0.08195359 0.22274308 0.6973931 0.79651105 0.3384423 0.5541741 0.60913694 0.11046143 0.7458725 0.1384661 0.7715413 0.34515342 0.301858 0.37151086 0.83048487 0.95463926 0.3078755 0.44282746 0.15716332 0.6104344 0.9658088 0.983976 0.9423979 0.94618285 1.001431 0.501199 0.9581948 0.1046057 0.9684777 0.766831 0.28662357 0.23676634 0.9742875 0.33788764 0.19063273 0.859799 0.47100848 0.37116095 0.11618819 0.31927577 0.13510212 0.98332024 0.24861254 0.9156055 0.9705818 0.93859595 0.25979528 0.6709827 0.6430121 0.96920294 0.98922735 0.1332146 0.9730394 0.9900741 0.9792404 0.6133692 0.9566109 0.38244438 0.6167695 0.6194038 0.4230758 0.5812086 0.98492885 0.95988905 0.96638143 0.6552822 0.95254976 0.95320314 0.33100057 0.14336959 0.15652771 0.20886818 0.13409692 0.9935843 0.77358496 0.65991676 0.26668385 0.98532397 0.9928422 0.5742656 0.13056873 0.8847005 0.1255912 0.21507831 0.9229842 0.10622708 0.8656443 0.20110266 0.88761157 0.94561636 0.3301112 0.98461807 0.23270014 0.5954498 0.96808904 0.17934097 0.98950607 0.9325035 0.97499657 0.8738652 0.6806155 0.95022196 0.98161185 0.3143995 0.11067381 0.96562856 0.96314186 0.8141581 0.97733825 0.92926216 0.99643946 0.9783236 0.44274423 0.9624729 0.5302121 0.9856205 0.42674586 0.8913692 0.5786059 0.40321144 0.23011269 0.8344016 0.8662287 0.3951003 0.12271074 0.9822601 0.98576945 0.08965831 0.23458704 0.8308076 0.5976968 0.9018738 0.64758146 0.9942369 0.98407173 0.88114756 0.98784196 0.71965104 0.9768777 0.97164696 0.8124566 0.18265249 0.07274331 0.1582355 0.93094766 0.8057281 0.16951917 0.96500504 0.11236409 0.20663357 0.18963169 0.9164737 0.13590461 0.12270089 0.37726036 0.96880966 0.20889229 0.95033556 0.6487601 0.41053978 0.9963639 0.63928187 0.8139373 0.12891752 0.9971641 0.95218897 0.33437562 0.06403207 0.63105154 0.6764621 0.5572445 0.0968667 0.9999528 0.99365973 0.8747377 0.8860685 0.9603575 0.76183 0.13849731 0.6251325 0.81368667 0.5709297 0.11879976 0.6601202 0.94263464 0.15936175 0.5716257 0.922511 0.24030536 0.9595926 0.5674498 0.29708254 0.3931772 0.12726173 0.9858892 0.9569603 0.97806525 0.9719149 0.9875281 0.9084232 0.23101103 0.94993997 0.9570631 0.4568754 0.19421495 0.67033285 0.11897636 0.12902626 0.9566067 0.35129988 0.9625631 0.6133013 0.8180542 0.29817158 0.99673927 0.99112725 0.74779296 0.299515 ] [1 0 1 1 0 0 1 1 0 0 0 1 1 1 1 1 0 0 0 1 0 1 1 0 1 1 0 0 1 0 1 1 0 0 1 0 0 1 0 1 0 1 0 0 1 0 1 1 0 0 1 0 0 0 1 0 1 1 0 0 0 0 0 0 1 0 1 1 1 1 1 0 1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 0 0 0 1 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 0 1 1 1 0 1 1 1 1 0 1 0 1 0 0 1 0 0 1 0 0 1 0 1 1 0 0 1 0 0 0 1 1 0 1 1 0 0 0 0 0 1 0 1 1 0 0 1 1 1 0 1 0 0 1 1 1 1 0 0 0 1 1 1 0 0 1 1 1 0 0 0 1 1 1 0 1 0 1 1 0 1 0 1 0 1 1 0 0 0 1 0 0 0 0 1 1 1 1 0 0 0 1 0 0 1 1 0 0 0 0 0 0 1 0 1 0 0 1 1 0 0 1 0 1 1 0 1 0 0 0 0 0 1 0 0 0 1] (100, 248)
print(all_y_pred)
[[0.2956851 0.1307303 0.8932625 ... 1.0380028 0.7612098 0.30560386] [0.21882373 0.2946883 1.0059178 ... 0.88635397 0.7102764 0.29325962] [0.41497844 0.31902927 0.9113966 ... 0.9928802 0.63424313 0.26908088] ... [0.2796316 0.22994864 0.96426916 ... 0.9784926 0.7533461 0.32739812] [0.24358472 0.27926666 0.97475976 ... 0.96031356 0.7478002 0.28509638] [0.30003503 0.25516906 0.95309937 ... 0.96751046 0.77303225 0.27185297]]
# Create a figure and axes
fig, ax = plt.subplots(1, 1, figsize=(14, 8))
# Plot predictions
ax.plot(np.arange(1, y_pred.shape[0] + 1), y_pred, label='Predictions', color='salmon', marker='.')
# Plot true values
ax.plot(np.arange(1, y.shape[0] + 1), y, label='True Values', color='lightseagreen', marker='.')
# Set limits, labels, grid, and legend
ax.set_ylim([0, 1])
ax.set_xlabel('Test Engine Units', fontsize=16)
ax.set_ylabel('RUL', fontsize=16)
ax.grid(True)
ax.legend()
plt.savefig(f'figures\\Responses_without_Uncertainty_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
# Show the plot
plt.show()
plt.plot(np.arange(1,len(T)+1), T, label= 'Train loss')
plt.plot(np.arange(1,len(V)+1), V, label = 'Validation loss')
plt.legend()
plt.grid()
plt.xlabel('Epochs')
plt.ylabel('MSE')
plt.savefig(f'figures\\Training_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show
<function matplotlib.pyplot.show(close=None, block=None)>
fig, ax = plt.subplots(figsize = (12,8))
def animate(i):
ax.clear()
line1, = ax.plot(np.arange(1,i+1), T[:i], label = 'train_loss')
line2, = ax.plot(np.arange(1,i+1), V[:i], label = 'val_loss')
ax.legend()
ax.grid(True)
ax.set_xlim(0,101)
ax.set_ylim(0,4000)
ax.set_xlabel('epochs')
ax.set_ylabel('MSE')
return line1, line2
animation = FA(fig, animate, np.arange(1,len(T)+1), interval = 50)
%time animation.save('animation3.gif', writer='imagemagick', fps=20)
plt.close(fig)
MovieWriter imagemagick unavailable; using Pillow instead.
CPU times: total: 21.4 s Wall time: 22 s
# Count total parameters
total_params = sum(p.numel() for p in model.parameters())
print(f"Total number of parameters: {total_params}")
# Count trainable parameters only
trainable_params = sum(p.numel() for p in model.parameters() if p.requires_grad)
print(f"Trainable parameters: {trainable_params}")
Total number of parameters: 22411 Trainable parameters: 22411
from torchinfo import summary
summary(model, input_size=(64, window, n_features))
=============================================================================================== Layer (type:depth-idx) Output Shape Param # =============================================================================================== RegressorModel [64, 1] -- ├─DEQFixedPoint: 1-1 [64, 60, 20] -- │ └─DELayer: 2-1 [64, 60, 20] -- │ │ └─Conv1d: 3-1 [64, 30, 20] 450 │ │ └─Conv1d: 3-2 [64, 30, 20] 2,700 │ │ └─Conv1d: 3-3 [64, 30, 20] 2,700 │ │ └─GroupNorm: 3-4 [64, 60, 20] 120 │ │ └─DualInputAttention: 3-5 [64, 60, 20] 1,680 │ │ └─GroupNorm: 3-6 [64, 60, 20] 120 │ └─DELayer: 2-2 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-7 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-8 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-9 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-10 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-11 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-12 [64, 60, 20] (recursive) │ └─DELayer: 2-3 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-13 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-14 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-15 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-16 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-17 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-18 [64, 60, 20] (recursive) │ └─DELayer: 2-4 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-19 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-20 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-21 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-22 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-23 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-24 [64, 60, 20] (recursive) │ └─DELayer: 2-5 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-25 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-26 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-27 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-28 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-29 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-30 [64, 60, 20] (recursive) │ └─DELayer: 2-6 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-31 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-32 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-33 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-34 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-35 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-36 [64, 60, 20] (recursive) │ └─DELayer: 2-7 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-37 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-38 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-39 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-40 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-41 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-42 [64, 60, 20] (recursive) │ └─DELayer: 2-8 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-43 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-44 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-45 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-46 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-47 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-48 [64, 60, 20] (recursive) │ └─DELayer: 2-9 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-49 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-50 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-51 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-52 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-53 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-54 [64, 60, 20] (recursive) │ └─DELayer: 2-10 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-55 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-56 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-57 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-58 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-59 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-60 [64, 60, 20] (recursive) │ └─DELayer: 2-11 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-61 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-62 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-63 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-64 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-65 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-66 [64, 60, 20] (recursive) │ └─DELayer: 2-12 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-67 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-68 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-69 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-70 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-71 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-72 [64, 60, 20] (recursive) │ └─DELayer: 2-13 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-73 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-74 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-75 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-76 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-77 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-78 [64, 60, 20] (recursive) │ └─DELayer: 2-14 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-79 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-80 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-81 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-82 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-83 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-84 [64, 60, 20] (recursive) │ └─DELayer: 2-15 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-85 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-86 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-87 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-88 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-89 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-90 [64, 60, 20] (recursive) │ └─DELayer: 2-16 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-91 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-92 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-93 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-94 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-95 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-96 [64, 60, 20] (recursive) │ └─DELayer: 2-17 [64, 60, 20] (recursive) │ │ └─Conv1d: 3-97 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-98 [64, 30, 20] (recursive) │ │ └─Conv1d: 3-99 [64, 30, 20] (recursive) │ │ └─GroupNorm: 3-100 [64, 60, 20] (recursive) │ │ └─DualInputAttention: 3-101 [64, 60, 20] (recursive) │ │ └─GroupNorm: 3-102 [64, 60, 20] (recursive) ├─AdaptiveAvgPool1d: 1-2 [64, 60, 1] -- ├─BNN: 1-3 [64, 1] -- │ └─Linear: 2-18 [64, 120] 7,320 │ └─Dropout: 2-19 [64, 120] -- │ └─Linear: 2-20 [64, 60] 7,260 │ └─Dropout: 2-21 [64, 60] -- │ └─Linear: 2-22 [64, 1] 61 =============================================================================================== Total params: 22,411 Trainable params: 22,411 Non-trainable params: 0 Total mult-adds (Units.MEGABYTES): 130.32 =============================================================================================== Input size (MB): 0.08 Forward/backward pass size (MB): 78.43 Params size (MB): 0.09 Estimated Total Size (MB): 78.60 ===============================================================================================
3. Prediction on Test Set¶
model.load_state_dict(torch.load(f'saved_models\\DEM_{file_name}', map_location=torch.device(device)))
C:\Users\spiro\AppData\Local\Temp\ipykernel_9692\4004060783.py:1: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature.
model.load_state_dict(torch.load(f'saved_models\\DEM_{file_name}', map_location=torch.device(device)))
<All keys matched successfully>
Phm08 Score Metric
def calculate_phm08_score(true_rul, predicted_rul, alpha=13, beta=10):
"""
Calculate the PHM08 score metric for RUL prediction.
Parameters:
true_rul (array-like): Array of true RUL values.
predicted_rul (array-like): Array of predicted RUL values.
alpha (float): Scaling factor for early predictions (default=13).
beta (float): Scaling factor for late predictions (default=10).
Returns:
float: Total PHM08 score.
"""
true_rul = np.array(true_rul)
predicted_rul = np.array(predicted_rul)
if true_rul.shape != predicted_rul.shape:
raise ValueError("Shape mismatch between true and predicted RUL arrays")
d = (predicted_rul - true_rul) * MAX_RUL
score = np.where(d < 0, np.exp(-d / alpha) - 1, np.exp(d / beta) - 1)
return np.sum(score)
num_samples = 100
def test(num_samples=num_samples):
model.train() # Enable stochastic behavior (dropout during inference)
total_loss_MSE = 0.0
total_loss_L1 = 0.0
total_ASUE = 0.0
total_phm08_score = 0.0
all_stochastic_predictions = []
all_pred_values = []
all_true_values = []
all_uncertainties = [] # Collect uncertainties
with torch.no_grad():
for X, y in testloader:
X, y = X.to(device).float(), y.to(device).float()
stochastic_predictions = []
# Perform stochastic forward passes
for _ in range(num_samples):
y_pred = model(X, training=True).squeeze()
stochastic_predictions.append(y_pred.cpu().numpy())
# Convert to NumPy array
stochastic_predictions = np.array(stochastic_predictions) # Shape: (num_samples, batch_size)
# Calculate mean and uncertainty with std
mean_pred = torch.tensor(np.mean(stochastic_predictions, axis=0)).to(device)
uncertainty_batch = np.std(stochastic_predictions, axis=0) # Variance as uncertainty
# Compute batch losses
loss_MSE = torch.mean((mean_pred - y) ** 2).item()
loss_L1_val = loss_L1(mean_pred, y).item()
ASUE = torch.mean(torch.relu(y - mean_pred)).item()
# Compute PHM08 score
phm08_score = calculate_phm08_score(y.cpu().numpy(), mean_pred.cpu().numpy())
# Aggregate results
total_loss_MSE += loss_MSE * X.size(0)
total_loss_L1 += loss_L1_val * X.size(0)
total_ASUE += ASUE * X.size(0)
total_phm08_score += phm08_score
# Collect predictions, true values, and uncertainties
all_stochastic_predictions.append(stochastic_predictions)
all_pred_values.append(mean_pred.cpu().numpy())
all_true_values.append(y.cpu().numpy())
all_uncertainties.append(uncertainty_batch)
# Average metrics
total_samples = len(testloader.dataset)
avg_loss_MSE = total_loss_MSE / total_samples
avg_loss_L1 = total_loss_L1 / total_samples
avg_ASUE = total_ASUE / total_samples
avg_phm08_score = total_phm08_score
# Concatenate all predictions, uncertainties, and true values
all_pred_values = np.concatenate([x.flatten() for x in all_pred_values], axis=0)
all_stochastic_predictions = np.concatenate(all_stochastic_predictions, axis=1)
all_true_values = np.concatenate(all_true_values, axis=0)
all_uncertainties = np.concatenate(all_uncertainties, axis=0) # Ensure consistency
return avg_loss_MSE, avg_loss_L1, avg_ASUE, avg_phm08_score, all_uncertainties, all_stochastic_predictions, all_pred_values, all_true_values
# Perform testing with Monte Carlo Dropout
mse, l1, asue, score, uncertainty, all_y_pred, y_pred, y = test()
# Calculate RMSE and scale by RUL max value
rmse = rul_max * np.sqrt(mse)
# Print results
print(f"Test Results:\n"
f"RMSE: {round(rmse, 2)}\n"
f"Score: {round(score, 2)}\n"
f"L1 Loss: {round(l1, 2)}\n"
f"ASUE: {round(asue, 2)}")
# Mean uncertainty
mean_uncertainty = np.mean(uncertainty)
print(f"Mean Uncertainty: {round(mean_uncertainty, 4)}")
# Optional: Visualize predictions, true values, and uncertainties
import matplotlib.pyplot as plt
plt.figure(figsize=(14, 8))
# True vs predicted RUL
plt.plot(range(len(y)), y, label="True Values", c="lightseagreen", marker=".", alpha=0.7)
plt.plot(range(len(y_pred)), y_pred, label="Predicted Values", c="salmon", marker=".", alpha=0.7)
# Uncertainty as error bars
plt.fill_between(
range(len(y_pred)),
y_pred - uncertainty,
y_pred + uncertainty,
color="salmon", alpha=0.2, label="Uncertainty"
)
# Add labels and legend
plt.xlabel("Sample Index", fontsize=14)
plt.ylabel("RUL", fontsize=14)
plt.title("True vs Predicted RUL with Uncertainty", fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.savefig(f'figures\\Responses_with_Uncertainty_for_dataset_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
Test Results: RMSE: 11.82 Score: 648.08 L1 Loss: 0.07 ASUE: 0.03 Mean Uncertainty: 0.05570000037550926
values = np.array([y, y_pred, uncertainty]).T
df_values = pd.DataFrame(values)
values = df_values.sort_values(0, ascending=False).to_numpy()
y, y_pred, uncertainty = np.array(values).T
# y, y_pred, uncertainty = MAX_RUL * np.array(values).T
# uncertainty = MAX_RUL * uncertainty
plt.hist(uncertainty, bins=50, color="blue", alpha=0.7, label="Uncertainty Distribution")
plt.xlabel("Uncertainty")
plt.ylabel("Frequency")
plt.title("Histogram of Prediction Uncertainty")
plt.legend()
plt.grid(True)
plt.show()
# confidence = 0.99 # Confidence level
# zeta = 2.576
# # Calculate confidence interval
# standard_error = uncertainty / np.sqrt(num_samples) # Standard error of the mean
# margin_of_error = zeta * standard_error # Margin of error
confidence = 0.95 # Confidence level
margin_of_error = confidence * uncertainty # Margin of error
plt.figure(figsize=(14, 8))
# True Values (more prominent line)
plt.plot(range(len(y)), y, label="True RUL", c="lightseagreen", marker="o", linewidth=2, alpha=0.8)
# Predicted Values (smaller, less opaque points)
plt.scatter(range(len(y_pred)), y_pred, label="Predicted RUL", c="salmon", edgecolor="black", s=30, alpha=0.7, zorder=3)
# Vertical Lines (improved clarity)
for i in range(len(y)):
plt.plot([i, i], [y_pred[i], y[i]], c="gray", linestyle="--", linewidth=0.8, alpha=0.7)
# Confidence Interval (adjusted opacity, more descriptive label)
plt.fill_between(
range(len(y_pred)),
y_pred - margin_of_error,
y_pred + margin_of_error,
color="salmon",
alpha=0.3, # Reduced opacity
label=f"{int(confidence * 100)}% Prediction Interval" # More accurate label
)
# Labels and Title (improved clarity)
plt.xlabel("Testing Sample Index", fontsize=14) # more context
plt.ylabel("Remaining Useful Life (cycles)", fontsize=14) # Full RUL
plt.title("Comparison of True and Predicted Remaining Useful Life (RUL)", fontsize=16) # more descriptive title
# Legend and Grid
plt.legend(fontsize=12)
plt.grid(True, alpha=0.5) # reduced grid alpha
# Save and Show
plt.savefig(f'figures/True_vs_PredictedRUL_with_Confidence_Interval_{file_name}.png', dpi=300, bbox_inches="tight")
plt.show()
import matplotlib as mpl
mpl.rcParams['font.family'] = 'serif' # Use serif font for better readability
mpl.rcParams['font.size'] = 12 # Consistent font size
mpl.rcParams['axes.linewidth'] = 1 # Thicker axes lines
mpl.rcParams['xtick.major.size'] = 5 # Ticks adjustments
mpl.rcParams['ytick.major.size'] = 5
mpl.rcParams['xtick.major.pad'] = 5
mpl.rcParams['ytick.major.pad'] = 5
mpl.rcParams['axes.labelsize'] = 14
mpl.rcParams['axes.titlesize'] = 16
mpl.rcParams['legend.fontsize'] = 12
mpl.rcParams['figure.dpi'] = 300 # High DPI for publication quality
mpl.rcParams['savefig.bbox'] = 'tight' # Tight bounding box
mpl.rcParams['axes.grid'] = True # Enable grid by default
mpl.rcParams['grid.alpha'] = 0.5 # Grid alpha
mpl.rcParams['lines.linewidth'] = 1.5 # Line width
mpl.rcParams['scatter.marker'] = 'o' # Marker
plt.figure(figsize=(8, 6)) # Adjust figure size for publication
# True Values (more prominent line)
plt.plot(range(len(y)), y, label="True RUL", c="lightseagreen", linewidth=2)
# Predicted Values (smaller scatter points)
plt.scatter(range(len(y_pred)), y_pred, label="Predicted RUL", c="salmon", edgecolor="black", alpha=0.7, zorder=3, s=15) # reduced scatter point size to 15
# Vertical Lines (improved clarity)
for i in range(len(y)):
plt.plot([i, i], [y_pred[i], y[i]], c="gray", linestyle="--", linewidth=0.8, alpha=0.7)
# Confidence Interval (adjusted opacity, more descriptive label)
plt.fill_between(
range(len(y_pred)),
y_pred - margin_of_error,
y_pred + margin_of_error,
color="salmon",
alpha=0.3, # Reduced opacity
label=f"{int(confidence * 100)}% Prediction Interval" # More accurate label
)
# Labels and Title (improved clarity)
dataset_name = file_name[:-4]
plt.xlabel("Testing Sample Index") # Removed unnecessary context from label
plt.ylabel("Remaining Useful Life (RUL)") # Removed unnecessary context from label
plt.title(f"Comparison of True and Predicted RUL for {dataset_name}") # Removed unnecessary context from title
# Legend and Grid
plt.legend()
# Save and Show
plt.savefig(f'figures/True_vs_PredictedRUL_{file_name}.png')
plt.show()
# Calculate errors and mean predictions
errors = np.abs(y_pred.mean(axis=0) - y)
# Scatter plot with enhancements
plt.figure(figsize=(8, 6))
sc = plt.scatter(uncertainty, errors, c=errors, cmap='viridis', alpha=0.7, edgecolor='k')
# Add colorbar
cbar = plt.colorbar(sc)
cbar.set_label('Error Magnitude')
# Add grid and labels
plt.grid(alpha=0.3)
plt.xlabel('Uncertainty', fontsize=12)
plt.ylabel('Error', fontsize=12)
plt.title('Uncertainty vs. Error', fontsize=14)
# Optional: Fit a trend line
z = np.polyfit(uncertainty, errors, 1)
p = np.poly1d(z)
plt.plot(uncertainty, p(uncertainty), color='red', linestyle='--', label='Trend Line')
# Add legend
plt.legend()
# Show the plot
plt.tight_layout()
plt.savefig(f'figures\\UncertaintyVSError_{file_name}.png', dpi=300, bbox_inches='tight')
plt.show()
# Compute statistics
correlation = np.corrcoef(uncertainty, errors)[0, 1]
print(f"Correlation between uncertainty and error: {correlation:.2f}")
# Evaluate performance in high/low uncertainty regions
threshold = np.percentile(uncertainty, 75) # Upper 25% uncertainty
high_uncertainty_mask = uncertainty > threshold
low_uncertainty_mask = uncertainty <= threshold
high_uncertainty_error = errors[high_uncertainty_mask].mean()
low_uncertainty_error = errors[low_uncertainty_mask].mean()
print(f"Mean error in high uncertainty regions: {high_uncertainty_error:.2f}")
print(f"Mean error in low uncertainty regions: {low_uncertainty_error:.2f}")
Correlation between uncertainty and error: 0.79 Mean error in high uncertainty regions: 0.49 Mean error in low uncertainty regions: 0.25
4. Prediction on Test Set Without Dropout in the Inference¶
model.load_state_dict(torch.load(f'saved_models\\DEM_{file_name}', map_location=torch.device(device)))
C:\Users\spiro\AppData\Local\Temp\ipykernel_9692\4004060783.py:1: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature.
model.load_state_dict(torch.load(f'saved_models\\DEM_{file_name}', map_location=torch.device(device)))
<All keys matched successfully>
def calculate_phm08_score(true_rul, predicted_rul, alpha=13, beta=10):
"""
Calculate the PHM08 score metric for RUL prediction.
Parameters:
true_rul (array-like): Array of true RUL values.
predicted_rul (array-like): Array of predicted RUL values.
alpha (float): Scaling factor for early predictions (default=13).
beta (float): Scaling factor for late predictions (default=10).
Returns:
float: Total PHM08 score.
"""
true_rul = np.array(true_rul)
predicted_rul = np.array(predicted_rul)
if true_rul.shape != predicted_rul.shape:
raise ValueError("Shape mismatch between true and predicted RUL arrays")
d = (predicted_rul - true_rul) * MAX_RUL
score = np.where(d < 0, np.exp(-d / alpha) - 1, np.exp(d / beta) - 1)
return np.sum(score)
def test_without_dropout():
model.eval() # Disable stochastic behavior (dropout during inference)
total_loss_MSE = 0.0
total_loss_L1 = 0.0
total_ASUE = 0.0
total_phm08_score = 0.0
all_pred_values = []
all_true_values = []
with torch.no_grad():
for X, y in testloader:
X, y = X.to(device).float(), y.to(device).float()
# Perform a single deterministic forward pass
y_pred = model(X).squeeze()
# Compute batch losses
loss_MSE = torch.mean((y_pred - y) ** 2).item()
loss_L1_val = loss_L1(y_pred, y).item()
ASUE = torch.mean(torch.relu(y - y_pred)).item()
# Compute PHM08 score
phm08_score = calculate_phm08_score(y.cpu().numpy(), y_pred.cpu().numpy())
# Aggregate results
total_loss_MSE += loss_MSE * X.size(0)
total_loss_L1 += loss_L1_val * X.size(0)
total_ASUE += ASUE * X.size(0)
total_phm08_score += phm08_score
# Collect predictions and true values
all_pred_values.append(y_pred.cpu().numpy())
all_true_values.append(y.cpu().numpy())
# Average metrics
total_samples = len(testloader.dataset)
avg_loss_MSE = total_loss_MSE / total_samples
avg_loss_L1 = total_loss_L1 / total_samples
avg_ASUE = total_ASUE / total_samples
avg_phm08_score = total_phm08_score
# Concatenate all predictions and true values
all_pred_values = np.concatenate([x.flatten() for x in all_pred_values], axis=0)
all_true_values = np.concatenate(all_true_values, axis=0)
return avg_loss_MSE, avg_loss_L1, avg_ASUE, avg_phm08_score, all_pred_values, all_true_values
# Perform testing without Monte Carlo Dropout
mse, l1, asue, score, y_pred, y = test_without_dropout()
# Calculate RMSE and scale by RUL max value
rmse = rul_max * np.sqrt(mse)
# Print results
print(f"Test Results:\n"
f"RMSE: {rmse:.2f}\n"
f"Score: {score:.2f}\n"
f"L1 Loss: {l1:.2f}\n"
f"ASUE: {asue:.3f}")
Test Results: RMSE: 11.81 Score: 648.47 L1 Loss: 0.07 ASUE: 0.029
5. Methods for Weighting Predictions by Uncertainty¶
# Perform testing with Monte Carlo Dropout
mse, l1, asue, score, uncertainty, all_y_pred, y_pred, y = test()
loss_MSE = np.mean((y_pred - y) ** 2).item()
rmse = rul_max * np.sqrt(mse)
print(f'RMSE:{rmse:.2f}')
RMSE:11.75
all_y_pred_mean = all_y_pred.mean(axis=0)
weights = abs(all_y_pred - all_y_pred_mean)
weights = np.exp(-weights)
# weights = 1 / weights
weights /= weights.sum(axis=0)
weighted_predictions = np.sum(weights * all_y_pred, axis=0)
mse = np.mean((weighted_predictions - y) ** 2).item()
rmse = rul_max * np.sqrt(mse)
print(f'RMSE:{rmse:.2f}')
RMSE:11.76
Threshold In Uncertainty
# Define a threshold based on deviations (e.g., 1.5x the mean deviation)
threshold = 1.5 * weights.mean(axis=0)
# Mask predictions with deviations exceeding the threshold
valid_predictions = np.where(weights <= threshold, all_y_pred, np.nan)
# Compute predictions as the mean of valid predictions
weighted_predictions = np.nanmean(valid_predictions, axis=0)
mse = np.mean((weighted_predictions - y) ** 2).item()
rmse = rul_max * np.sqrt(mse)
print(f'RMSE:{rmse:.2f}')
RMSE:11.75